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INTRODUCTION 


It  is  the  responsibility  of  the  practitioner,  that  is  the 
engineer  or  scientist,  to  develop  the  simplest  model  that 
explains  or  predicts  variables  of  interest.  The  implications  of 
not  enough  model  detail  are  much  better  appreciated  than  the 
implications  of  too  much  model  detail.  If  the  entire  system  model 
identification  and  parameter  estimation  proceed  interactively 
with  experimentation  and  data  acquisition,  the  appropriate  model 
should  result. 

When  this  does  not  happen,  the  major  likely  reason  is:  the 
model  of  the  overall  system  is  typically  formed  from  many  small 
sub-models  and  parameters  are  estimated  from  data  on  the 
submodels.  Consider,  for  example,  the  kinetic  model  for 
combustion  of  a  hydrocarbon.  Hundreds  of  reactions  may  be 
theoretically  possible  to  occur  although  the  vast  majority  may 
not  need  to  be  considered  in  order  to  predict,  say,  the  overall 
burning  rate.  However,  typically  the  big  model  is  formed  and  then 
data  taken  from  experiments  on  many  small  submodels,  here 
consisting  of  the  individual  reactions.  Data  may  be  retrievable 
for  this  single  reaction  or  for  a  small  sub-set  of  reactions,  but 
it  may  not  be  important  or  desireable  to  do  so. 

The  numerical  penalty  for  too  much  model  detail  is  too  much 
computation  time,  either  because  the  system  is  larger  than  needed 
or  because  the  system  is  stiff,  or  both.  I  mean  here  to  use  this 


term  stiff  as  the  Webster  anonym  "subborn";  thus  included  in  this 
review  are  highly  oscillatory  systems,  "stiff"  two-point 
boundary  value  problems,  discontinuities,  and  stiff  boundary 
vaule  problems.  For  a  precise  definition  of  stiffness  see 
Shampine( 1982) ;  also  the  versions  of  Watson( 1976) , 
Robertson( 1976 ) ;  Scherer( 1976) ;  Lambert ( 1980) ;  Seider ( 1980) ;  and 
Gear(1980).  Guderley( 1975)  defines  stiff  two-point  boundary  value 
problems. 

There  is  another  fundamental  problem  with  the  identification 
and  estimation  staqe  beyond  the  inappropriate  model  detail  that 
results:  the  basis  for  the  derived  parameters  is  wrong.  Pretend 
that  you  must  design  a  process  for  separating  two  liquids  by 
boiling  them,  condensing  the  resulting  vapor  and  then  boiling 
again  in  stages  (distillation).  A  major  parameter  need  for  your 
design  equations  would  be  the  relative  volatilities  of  the 
components  as  a  function  of  composition.  This  can  be  done  in  a 
simple  bomb  experiment  in  the  laboratory  followed  by  least- 
squares  fitting  the  data  to  an  assumed  function;  but  in  doing 
this  you  have  obtained  results  optimal  at  best  to  your  rather 
arbitary  least  squares  criterion.  What  you  really  want  is  to 
minimize  the  uncertainty  in  the  design  of  your  process,  and  this 
should  be  expressed  at  the  parameter  estimation  stage. 

Let  us  consider  now  the  case  where  we  in  fact  do  a  "systems" 
identification  and  parameter  estimation  on  a  detailed  dynamic 


model.  We  then  face  the  inverse  problem  of  stiff  numerical 
solution  -  stiff  parameter  estimation.  This  not  only  can  present 
huge  experimental  requirements  but  the  need  also  for  special 
mathematical  and  numerical  considerations  (  Aiken,  1982;  Aiken 
and  Venkateschwaran,  1982). 

It  is  implicitly  assumed  in  most  all  the  papers  in  this 
review  that  a  stiff  model  is  desirable  from  a  predictive 
standpoint  and  that  the  parameters  are  appropriate  and  correct 
for  the  model.  If  this  is  not  true,  then  model  simplification  may 
be  in  order.  The  steady-state  approximation  applied  (properly)  a 
posteriori,  or  after  the  model  and  data,  is  a  mathematical 
approximation (Aiken  and  Lapidus,  1974;  Aiken  and  Lapidus,  1975); 
but  the  same  reasoning  applied  a  prori  is  just  a  modelling 
decision.  So  this  approximation  can  be  a  very  reasonable,  and 
powerful,  strategic  approach  before  attacking  with  the  stiff 
artillery. 

There  is  considerable  overlap  amonst  the  stiff  application 
areas  listed  below.  An  effort  is  made  in  each  area  to  explain  why 
the  models  are  stiff  and  what  some  of  the  researchers  did  to 
overcome  the  numerical  simulation  difficulties.  Industrial 
laboratories  are  explicitly  noted. 


ATMOSPHERIC 


Description  of  atmospheric  phenomena  involves  transport  with 
chemical  reaction;  stiffness  can  occur  because  the  time  scales  of 
the  reactions  are  much  smaller  than  times  for  transportation  over 
distances  of  interest.  The  transport  equations  themselves  can  be 
stiff  because  these  distances  can  be  very  lonq  indeed;  and  the 
chemical  kinetic  rate  equations  are  usually  large  and  usually 
stiff.  Pollutant  formation  models  naturally  are  stiff  because 
highly  reactive  free  radical  transients  are  included  in  the  model 
as  they  are  intimately  related  to  formation  of  trace  quantities 
of  polutants. 

Several  very  effective  packages  have  been  developed  at 
Lawrence  Livermore  Rational  Laboratory  for  the  solution  of 
differential  equations  modelling  atmospheric  phenomena  (Byrne, 
1981;  Hindmarsh,  1982).  EPISODE  was  written  for  stiff  chemical 
kinetics  including  minor  species  in  the  upper  atmosphere.  It  has 
solved  successfully  such  problems  with  diurnally  varying  reaction 
rates.  EPISODE  has  been  modified  (EPISODEB)  for  problems  that 
also  include  transport.  The  modification  is  to  recognize  the 
banded  structure  of  the  large  Jacobian  resulting  from  the  method 
of  lines  and  thus  to  handle  this  matrix  more  efficiently.  If 
finite  elements  or  collocation-B>splines  are  used  on  the  partial 
differential  equations,  the  resultinq  form 

A(y,t)dy/dt«f (y,t) 


where  A  is  of  a  banded  form  is  best  handled  by  another  version  of 
EPISODE,  EPISODEIB.  Systems  of  this  type  include  differential- 
alqebraic  sets  if  A  has  one  or  more  zero  columns.  Carmichael  et 
al.(1980),  for  example,  applied  the  Galerkin  finite  element 
method  to  the  movement  of  pollutants  in  the  atmosphere. 

Brown(1980),  however,  used  EPISODE  for  solution  of  a  diurnal 
kinetics  example  and  judges  it  too  inefficient.  Here  some 
concentrations  are  small  during  the  night,  suddenly  increase  by 
orders  of  magnitude  with  the  first  coming  of  light,  vary  slowly 
during  the  day  with  maximum  values  around  noon,  and  then  drop 
sharply  at  sunset.  With  the  cycle  time  well-know.  Brown 
transforms  the  time  coordinate  to  stretch  it  during  the  times  of 
rapid  change;  in  the  transformed  time  the  step  size  is  more 
constant. 

Miller  et  al.(T978)  present  a  one-dimensional  model  of 
atmospheric  fluorocarbon-ozone  photochemistry  with  transport. 
Chang  et  al.(1974)  solve  over  14,000  ordinary  differential 
equations  to  describe  the  effect  of  the  SST  on  the  ozone  layer. 

Bottenheim  and  Strausz( 1980)  modeled  gas-phase  chemistry  in 
clean  air  as  a  prelude  to  including  polutants.  Gelinas  and 
Skewes-Cox( 1977)  explore  tropospheric  photochemical  mechanisms. 
Baldwin  et  al.(1977)  study  rate  parameter  estimates  in  a 
photochemical  smog  kinetics  model.  Atkinson  et  al.(1980)  did  a 
modeling  study  fo  the  gas  -phase  NOX-air  photoxidation  of  toluene 


and  the  cresols.  Pitts  and  Finlayson( 1975)  propose  various 
mechanisms  of  photochemical  air  pollution.  Preussner  and 
Brand(1981)  apply  a  semi-implicit  Euler  method  to  photochemical 
smoq  kinetics.  Kuhlman  et  al.<1978)  study  the  effect  of  CO  on 
sulfate  aersol  formation.  Wallace  et  al.(1980)  model 
photochemical  ozone  and  NO  formation.  Watkins( 1981 )  solves  an 
ionospheric  model  with  several  unknown  initial  conditions;  these 
are  chosen  so  to  avoid  the  initial  steep  transient (see  also 
Aiken,  1974,  for  problem  approximations  like  this). 
Isaacson( 1981 )  predicts  extremely  high  wind  speeds  at  ground 
level  on  the  downstream  side  of  a  mountain  range. 

Difficulties  with  the  use  of  steady-state  approximation  on 
the  st if f  reaction  rate  equations  has  been  pointed-out  by  Farrow 
and  Edelson( 1974) .  Reasons  could  include:  the  necessity  of 
including  in  the  model  radicals  which  very  directly  effect  the 
trace  variables  of  interest;  the  re-occuring  nature  of  the 
stiffness  on  a  daily  basis;  the  typically  large  size  of  the 
equation  set  makes  choice  of  stiff  variables  more  difficult. 
Dickinson  and  Gelinas( 1976)  performed  sensitivity  analysis  on 
these  types  of  equation  to  better  understand  what  reactions  are 
important  in  their  example.  Farrow  and  GraedeK 1977)  found  the 
steady-state  approximations  applied  to  some  species  but  not  to 
others  that  they  would  have  expected  to  work;  some  species  rates 
can  have  this  approximation  hold  at  various  diurnal  times. 


BIOLOGICAL 


Most  chemical  reactions  occuring  in  living  species  involve  a 
catalyst  that  helps  with  proper  orientation  of  the  big  floppy 
biological  molecules.  The  catalyst  is  called  an  enzyme.  These 
reactions  are  stiff  for  one  of  two  reasons:  1 .  There  is  typically 
a  large  scale  difference  between  the  concentration  of  the 
reactant,  called  the  substrate,  and  the  enzyme,  or  2.  There  is  a 
large  difference  in  the  rate  constants.  Aiken (1982)  explores  the 
validity  and  the  implications  of  the  so-called  Michaelis- 
Menten(M-M)  approximation  to  relieve  the  stiffness  and  presents  a 
number  of  other  approximations.  These  could  be  more  valid 
depending  upon  the  relative  magnitudes  of  the  rate  constants. 

It  has  been  historically  rather  rare  to  find  enzyme  models 
where  the  M-M  approximation  has  not  been  used.  The  reason  for 
this  seems  to  be  primarily  an  experimental  limitation:  the 
enzyme-substrate  complex,  typically  the  stiff  intermediate 
variable,  is  like  the  enzyme  itself  present  in  very  small 
quantities  and  may  not  differ  greatly  from  measureable 
characteristics  of  the  enzyme  or  the  substrate.  It  was  measured 
spectroscopically  at  about  mid-twentieth  century  and  allowed  the 
first  complete  model  of  enzyme  kinetics  with  all  parameters 
specified.  There  seems  to  be  a  growing  interest  in  obtaining  this 
(stiff)  model  detail  today  (Hiromi,  1982;  Kondo,  et  al.,  1980). 
Kinderlehrer  and  Ainsworth ( 1976)  have  written  a  program  to 


oredict  detailed  enzyme  mechanisms  with  all  intermediates. 

In  addition  to  enzyme  kinetics  stiffness  appears  in 
pharmocological  compartment  models  of  drug  response.  Gehring  and 
Blau(1977),  of  Dow  Chemical,  have  modeled  dose  response  to 
suspected  carcinogens  and  noted  the  initial  large  transients. 
Bloch  et  al.(1980)  found  the  processes  of  reversible  binding  of 
drug  to  protein  occurs  rapidly  on  the  time  scale  of  the  solution 
to  the  drug  disposition  in  the  body.  Jackson( 1980)  used  a  version 
of  GEAR  for  the  kinetic  simulation  of  anticancer  drug 
interactions.  Perelson  and  DeLisi(1980)  encountered  stiffness 
with  receptor  clustering  on  a  cell  surface.  Perelson( 1979)  found 
that  an  antigen  will  come  on  and  off  the  surface  of  a  cell  many 
times  before  cross-linking  occurs.  Karba  et  al.(1980)  use  hybrid 
computer  simulation  to  overcome  stiffness  in  drug  pharmoco 
kinetics. 

Loomis  et  al.(1979)  , in  modelling  crop  physiology,  found 
that  while  he  was  interested  in  crop  growth  over  a  period  of  many 
days,  he  had  to  use  a  time  scale  of  hours  to  include  diurnal 
events  or  a  scale  of  minutes  or  seconds  if  cellular  process  were 
to  be  included.  Chu  and  Berman(1974)  developed  a  program  using  an 
exponential  method  for  modelling  and  simulation  of  complex 
bioloqical  systems.  Gottwald  and  Wanner(1982)  compare  various 
methods  for  stiff  differential  equations  occuring  in  biology. 
Hunding( 1980)  came  across  stiffness  and  "chemical  histeresis"  in 


simulating  a  biological  cell  or  early  blastula. 

COMBUSTION 

This  is  a  special  area  of  reaction  with  transport  and  of 
chemical  kinetics.  Stiffness  can  come  from  the  time  scale 
differences  between  reaction  and  transport  or  from  the  stiff 
kinetic  system  or  from  both. 

Combustion  kinetics  could  be  defined  as  simply  describing 
oxidative-type  reactions,  but  the  primary  reactions  of  interest 
are  initiated  by  combination  of  oxygen  with  hydrocarbons.  Prom 
this  class  the  most  important  reactions  are  the  chain  branching 
type  or  autocatalytic  reactions  that  proceed  ever  faster  as  fuel 
is  consumed  -  explosions.  Most  combustion  applications  involve 
explosive  reactions  confined  in  space;  the  spacial  location  of 
most  intense  reactivity  is  termed  the  flame  heart.  Aiken(1982) 
explores  the  definition  of  the  term  "explosion”  and  develops 
criteria  for  the  degree  of  explosive  activity  for  the  oxygen- 
hydrogen  problem. 

.Stiffness  in  explosive  reactions  is  different  from  that  in 
most  other  kinetic  systems  in  that  fast  transients  do  not  occur 
only  for  initial  times  but  usually  appear  later  in  the  transient. 
Thus  codes  that  are  taking  giant  time  integation  steps  can  jump 
over  the  explosive  zone  where  most  of  the  action  is.  Another 
difference  is  that  the  stiff  variables  change  during  the 
transient.  Typically  at  least  three  major  time  scales  are 


important  in  these  systems:  initiation,  "explosion”  and 
termination.  For  all  these  reasons  stiffness  in  combustion 
kinetics  can  be  particularly  severe;  only  recently  has  the 
steady-state  approximation  been  successfully  and  automatically 
applied  to  these  reactions  (Aiken,  1982).  Hoppensteadt  et 
al.(1981)  presents  a  projection  method  that  focuses  on  the 
positive  eigenvalue  during  the  explosion.  Pratt(1979)  has 
investigated  special  methods  that  make  use  of  the  special 
structure  of  the  kinetic  equations.  See  also  White  and 
Seider( 1981 ) . 

Addition  of  the  spacial  variable  further  complicates  the 
potential  numerical  problems.  The  flow  in  which  combustion  occurs 
is  usually  turbulent  with  intermittent  changes  on  time  scales 
different  from  the  kinetic  time  scales.  Stochastic  methods  are 
promising  for  overcoming  this  stiffness  (Chorin,  1980).  Kansa 
(1981)  ,of  the  Lawrence  Livermore  National  Laboratory,  has 
combined  some  aspects  of  block  implicit  PDE  schemes  with  stiff 
ODE  solvers  for  this  problem;  this  was  partially  motivated  by  the 
experience  that  block  implicit  methods,  although  unconditionally 
stable  for  linear  equations,  were  found  to  have  inadequate 
stability  properties  for  the  severe  nonlinear  Arrhenius 
temperature  dependence  of  the  reaction  rate. 

Operator  splitting  techniques  are  often  used  for  the 


numerical  solution  of  multicomponent  gas  mixtures  undergoing 

r 

rapid  reaction.  They  offer  better  storage  economy  than  fully 
implicit  schemes  and  better  stability  properties  than  explicit 
schemes;  their  low  order  accuracy  is  appropriate#  and  they  allow 
flexibility  and  modularity  in  the  overall  numerical  attack. 
Karasalo  and  Kurylo  (1981)#  also  of  Lawrence  Livermore#  discuss 
in  this  context  ways  to  reduce  overall  computation  time  by 
improving  the  efficiency  of  the  stiff  kinetics  step  at  each  grid 
point  where  by  far  the  most  time  is  spent.  They  employ  a  stiff 
ODE  package  (like  GEAR)  but  with  the  following  modifications: 
Pirst,  they  avoid  repeated  evaluations  of  the  Jacobian  at  each 
mesh  point  by  storing  selected  historical  values  from  adjacent 
mesh  points.  Second#  they  allow  step  sizes  and  order  to  vary  more 
frequently  since  this  is  not  the  expensive  feature.  Thirdly#  they 
artifically  impose  that  concentrations  remain  non-negative  during 
prediction  and  corrector  iterations.  This  last  point  is  an 
important  one  as  it  has  been  found  that  small  negative 
concentrations  that  can  result  from  rounding  or  truncation  causes 
stability  probems  that  require  a  local  error  bound  smaller  than 
that  otherwise  sufficient  for  a  given  global  accuracy  request. 

Otey(1978)#  of  Sandia#  formulated  a  test  problem  with 
oombustive  reaction  and  diffusion  to  make  comparisons  among 
solutions  by  the  method  of  lines#  linearized  block  tridiagonal 


procedures,  and  other  techniques.  Re  found  the  block  diagonal 
implicit  procedure  was  by  far  the  best  in  all  variants  of  the 
test  problem,  including  ones  for  which  the  system  was  quite  stiff 
through  stiffness  in  the  kinetic  equations.  Wendt  et  al.(1979) 
and  Wendt (1982)  present  methods  for  solving  stiff  boundary-value 
problems  with  combustion  and  diffusion  in  a  pore.  The  reaction 
occurs  on  the  surface  of  the  pore.  A  variable  grid  mesh  is  used 
so  as  to  be  finer  at  the  pore  mouth  where  concentrations  change 
fastest.  Sundaresan  and  AmundsonC 1980)  also  report  the  very  stiff 
nature  of  this  problem. 

Oran  and  Boris(1981)  present  a  comprehensive  exaimination  of 
modelling  and  simulation  of  combustive  flow  problems.  They  also 
mention  that  kinetic  ODE  stiffness  within  a  spacial  model  cannot 
be  handled  adequately  for  realistic  models  because  the  usual 
matrix  inversions  would  be  too  expensive  for  a  large  number  of 
chemical  species  and  many  grid  points;  storage  is  also  an  obvious 
biq  problem.  These  authors  review  the  very  complex  nature  of 
modelling  and  simulating  turbulence  for  reactive  environments. 

Sandusky  et  al.(1979)  advocate  finite  element  techniques  for 
combustive  transport.  Westbrock( 1978)  offers  an  improvement  in  th 
operator  splitting  method.  George  and  Harris (1977)  lament 
stiffness  from  a  model  of  in  situ  oil  shale  retorting.  See  also 
Scaccia  and  Kennedy ( 1 974 )  ,  McDonald ( 1979) ,  Choi  and 
Churchill ( 1979) ,  and  Lowe  et  al.(1977). 


CONTROL 


Many  engineering  control  systems  can  be  modelled  by  the 

forms 

dx/dt  »  Ax  +  Bz 
e  dz/dt  »  Cz  +  Bu 

where  the  first  equation  represents  a  large  linear  multi-variable 
plant,  and  the  second  equation  represents  a  multivariable 
actuator,  e  is  a  small  parameter  that  indicates  a  fast  contoller 
response  in  comparison  to  plant  variable  time  scale.  The 
controller  might  be  electrical  and  the  plant  mechanical  or 
involving  fluid  transport. In  general  e  cannot  be  neglected  (set 
equal  to  zero)  since  the  presence  of  the  controller  dynamics  can 
govern  the  inherent  stability  of  the  plant-controller  complex; 
see  Porter ( 1976 ) ,  however,  for  some  conditions  for  which  this 
simplification  is  permissible.  Shimizu  et  al.(1980)  describes 
some  stiff  nonlinear  control  problems. 

This  singularly  perturbed  simple  linear  form  has  been 
exploited  by  Khalil  and  Kokotovic( 1980)  in  the  design  of 
decentralized  feedback  controlers.  Anderson ( 1980)  offers  a  time- 
varying  transformation  to  separate  the  fast  from  the  slow  modes. 
Dontchev( 1974)  explores  by  sensitivity  analysis  optimal  control 
systems  with  changes  in  order.  Womble(1976)  looks  at  some  further 
approximations  to  Ricatti  equations  having  fast  and  slow  modes. 

Any  differential  equation  can  be  considered  stiff  if  a 


solution  is  required  in  a  short  enough  computing  time.  Real  time 
aerospace  control  applications  can  have  this  feature.  Bulirsch 
and  8ranca(1974)  mention  for  an  optimal  deceleration  maneuver  an 
Apollo-type  vechicle  would  need  information  in  about  one  second 
and  would  allow  up  to  1%  error .Gear  (1977)  discusses  the  conflict 
between  real-time  and  software;  he  mentions  that  real-time 
implies  that  implicit  methods  cannot  be  used  in  the  usual  sense 
and  presents  some  of  the  semi-implicit  methods.  See  Hiestand  and 
George(1976)  for  other  stiff  aerospace  applications. 

Ojika  et  al.  (1979)  give  a  "time  decompostion"  algorithm  for 
a  stiff  two  point  boundary  value  problem  applied  to  nonlinear 
optimal  control  problems. 

DISPERSED  PHASES 

Consider  a  liquid  or  solid  phase  dispersed  as  droplets  or 
particles  in  another  gas,  liquid  or  solid  phase.  If  the  size 
distribution  of  the  dispersed  phase  is  broad,  stiffness  can 
result  from  models  that  include  heat  or  mass  transfer, 
particularly  complex  when  reaction  is  also  occuring  in  the 
dispersed  phase.  This  is  because  it  is  much  easier  to  transport 
mass  and  heat  to  smaller  sizes,  roughly  proportional  to  the 
reciprocal  of  the  effective  diameter.  For  example,  Kayihan( 1980) 
describes  a  model  and  solution  for  coal  devolatilization  in  which 
heat  is  transfered  preferentially  to  the  smaller  particles  that 


therefore  devolatize  relatively  rapidly  and  cause  severe 
numerical  simulation  problems.  As  usual ,  a  modelling  choice  must 
be  made  as  how  small  a  size  to  include.  The  size  distribution 
functions  typically  have  a  skewed  maximum  with  long  tails  - 
particularly  in  the  direction  of  the  smaller  sizes.  The  modeller 
must  often  determine  a  cut-off  point  for  small  particles;  smaller 
than  that  would  cause  numerical  problems  but  ideally  not 
contribute  significantly  to  the  solution  variables  of  interest. 

Wall  and  Anlansson( 1980)  use  a  version  of  the  GEAR  package 
to  solve  a  model  of  stepwise  micelle  association.  Lahey  et 
al.(1980)  also  use  this  package  for  modelling  bubbles  flowing 
through  a  nozzle. 

See  also  the  sections  in  this  review  on  heat  trar^se-  and 
reactors.  Bubble  columns  and  spray  reactors  can  exhibit  this  type 
of  stiffness  as  can  processing  crushed  shale  or  coal.  This 

problem  can  be  particularly  severe  in  the  case  of  in  situ 
solids  processing  for  which  the  particle  sizing  from  underground 
explosions,  for  example,  is  widely  varying. 

ELECTRONICS 

The  time  domain  analysis  of  electronic  circuits  requires  the 
solution  of  nonlinear  algebraic-differential  equations.  Implicit 
integration  methods  and  sparse  matrix  techniques  made  possible 
analysis  of  circuits  containing  hundreds  of  active 


devices 


Advances  in  large-scale  integrated  circuits  have  indicated 
potential  for  analysis  of  thousands  of  active  elements.  Hybrid 
method  simulation  is  an  interesting  concept  that  applies 
different  methods  to  sections  of  the  circuit  that  require 
different  accuracy,  but  effects  of  interaction  among  the 
subsystems  can  be  difficult  to  assess  a  priori.  The  concept  of 
"latency",  rather  like  a  temporary  steady  state  approximation, 
has  also  received  attention  in  this  area.  The  relationship 
between  latency  and  the  numerical  method  has  been  explored  in 
Rabbat  et  al.(1979),  of  IBM  Data  Systems  Division.  See  also 
Tadeusiewicz( 1981 ) . 

Power  system  dynamic  response  involves  the  solution  of  large 
dif ferential-alqebraic  equations.  The  differential  equations 
model  the  dynamics  of  the  machines  and  their  control  systems 
while  the  algebraic  equations  model  the  network  steady-state 
relationships.  Gross  and  Bergen{1977)  pursue  this  combination  by 
partitioning  the  set  into  a  non-stiff  part  and  a  stiff  part  with 
a  sparse  Jacobian  matrix. 

Resonant  circuits,  time  variant  and  time  invariant,  have 
been  studied  by  Ruehli  et  al.(1980),  of  IBM  T.J.  Watson  Research 
Center.  "A-contractive  arc"  methods  were  shown  to  perform  well 
for  both  types  of  circuits.  Oscillatory  nonlinear  circuits  with  a 
finite  number  of  continuous  derivatives  has  been  the  subject  of 


work  by  Hajj  and  Skelboe{ 1979) .  Zein(1980),  of  IBM  Data  Systems 
Division,  discusses  the  use  of  "sensitivity  circuits”  for  the 
transient  analysis  of  periodic  circuit  behavior. 

A  range  of  integration  algorthms  have  been  tested  on  some 
model  problems  of  large  ODE  sets  for  power  system  dynamics  by 
Humpage  et  al.(1974).  Methods  for  fast  contingency  analysis  at  a 
power  systems  control  center  are  advanced  by  Chamorro  et 
al.  ( 1981 ) . 

Alvarado ( 1979)  reports  some  results  on  stiff  transient 
stability  analysis  of  circuits.  Weaver  et  al.(1977)  give  a  stiff 
model  for  radiation-induced  bulk  electrical  conductivity  in 
insulators.  Covello  and  White{1977),  of  the  U.S.  Air  Force 
Weapons  Laboratory,  discuss  stiffness  when  investigating 
radiation  response  characteristics  of  networks. Charge  transfer  in 
a  nonlinear  stiff  model  of  charge-coupled  devices  has  been 
simulated  by  McKenna  and  Schryer( 1975) ,  from  Bell  Labs.  Warner 
and  Wilson( 1980) ,  also  working  at  Bell  Labs,  use  some  analytic 
transformations  to  help  lessen  the  stiffness  from  their  equations 
related  to  the  fabrication  of  narrow-channel  MOS  transistors. 
Gambart-Ducros  and  Maral(1980)  concern  themselves  with  the  stiff 
differential  equations  that  arise  from  some  computer  aided  design 
techniques;  see  also  Dietze  and  Reibiger( 1978 ) .  Von  Pragenau 
(1981)  reports  his  own  method  for  greatly  reducing  the 
computation  time  of  stiff  digitial  filters.  Stiff  nonlinear 


switching  circuits  has  been  attacked  by  Boness(1979)  and 
specifically  switching  serges  byTripathy  and  Rao(1978). 

See  the  section  in  this  review  on  control  systems  for  some 
further  stiff  applications  involving  electronics. 

FLUIDS 

Stiffness  occurs  in  spacial  coordinates  within  homogeneous 
fluids  with  sharp  changes  in  physical  properties  or  from  abrupt 
obstructions  in  the  flow  path.  Compressable  flow  with  compression 
and  rarefaction  (shock)  waves,  reflection,  flow  reversal  and 
choked  flow  all  can  lead  to  numerical  problems. 

The  method  of  characteristics  is  "characteristically"  used 
on  problems  of  inviscid  flow  because  it  naturally  handles 
discontinuous  derivatives  as  it  follows  waves  but  cannot  be  used 
on  viscous  shock  layer  equations.  Srivastava  et  al.(1979)  present 
a  finite  differencing  scheme  for  viscous  flow  past  blunted  cones, 
where  derivative  discontinuities  are  encountered  at  the  sphere- 
cone  juncture  point.  To  avoid  large  truncation  errors  associated 
with  these  points,  differencing  across  the  discontinuity  is 
carefully  avoided.  The  method  of  characteristics  can  become  too 
expensive  on  inviscid  flow  problems  to  follow  long  term 
transients  involving  shock  waves;  Carver(1980)  gives  a  spacial 
discretization  which  utilizes  the  directional  aspects  of  the 
method  of  characteristics. 


Rlottner( 1980) ,  from  Sandia  Laboratory,  used  a  variable  grid 
approach  to  solve  turbulent  boundary-layer  flows  that  involve 
jumps  in  viscosity  and  density.  Blottner  has  used  both  coordinate 
stretching  with  uniform  grid  in  the  stretched  coordinate  as  well 
as  discontinuous  grid  spacing  that  is  effective  for  discontinous 
changes  in  variables.  MacCormack  and  Paullay( 1974)  ,of  NASA  Ames, 
provide  a  study  on  the  effect  of  the  mesh  spacing  on  inviscid 
supersonic  shock  flows.  Stewart ( 1979) ,  of  the  Atomic  Energy 
Laboratory  in  France,  examines  a  model-oriented  numerical  method 
for  solving  flow  with  sharp  changes  in  phase  as  occurs  in  cooling 
water  superheated  locally.  Other  variable  mesh  approaches  are 
shown  by  the  Russian  group  of  Yanenko( 1979)  for  boundary  layer 
shear  flows. 

The  POEs  that  model  unsteady  flow  in  one,  two  or  three 
dimensions  can  yield  stiff  ODEs  when  spacially  discretized  -  the 
method  of  lines.  This  will  be  the  case  where  there  are  many 
spacial  grid  points  compared  to  the  time  steps  that  one  would 
like  to  use.  A  large  number  of  spacical  mesh  points  can  result 
from  tight  spacial  coupling  or  simply  from  "long"  distances  to  be 
covered.  The  method  of  lines  can  be  attractive  because  of  its 
programming  simplicity  but  is  not  as  efficient  as  finite 
differencing  for  problems  without  tight  spacial  coupling  (Kurtz 
et  al. ,  1978) . 

Nadala  and  Piacsek( 1977)  ,of  the  U.S.  Naval  Laboratory,  have 


studied  numerically  the  responce  of  oceans  to  weather  changes. 
They  avoid  small  time  steps  associated  with  fast  moving  surface 
gravity  waves  by  dividing  the  flow  into  baroclinic  and  barotropic 
vertically  averaged  modes;  the  baroclinic  waves  are  treated 
explicitly  and  the  barotropic  waves  implicitly  (still  computation 
times  reach  over  60  hours  on  a  large  computer). 

See  also  Coleman  et  al.(1977)  for  averaging  methods  applied 
to  stiff  circulatory  flow  and  Gersting( 1980)  for  the  Orr- 
Somerfield  flow  approached  as  an  initial  value  problem. 
Issacson( 1981 )  looks  at  the  mountain  wind  problem  and  suggests  a 
filtering  scheme  and  a  hybrid  method  for  handling  shocks  in  the 
atmosphere. 

Refer  to  sections  in  this  review  on  combustion ,  reactors , 
atmospheric  problems  and  general  reaction  -  diffusion  for  flow 
problems  with  reaction  coupled  to  transport. 

HEAT 

Stiffness  in  heat  transfer  originates  in  one  of  two  ways: 
sharp  changes  in  temperature  environment  or  large  differences  in 
the  rates  which  components  of  the  system  can  transfer  heat.  The 
first  problem  could  be  a  boundary-value  problem  with  the  sharp 
changes  represented  in  the  boundary  conditions.  A  realistic  model 
for  these  boundary  conditions  is  a  tough  problem  in  itself  since 
discontinuities  would  not  exist  in  nature;  step  changes  in 


temperature  can  result  in  infinite  heat  fluxes,  certainly  not 
observed.  The  other  class  of  stiff  heat  transfer  problems  arise 
from  differences  in  heat  capacities  or  by  size  differnces  among 
the  components.  For  this  category  see  also  the  section  on 
dispersed  phase  transport  where,  for  example,  heating  of  coal 
particles  for  pyrolysis  is  discussed. 

Krishnan  and  Sastri(1978)  solved  the  thermal  entry  length 
« 

problem  for  high  Prandtl  numbers,  that  is,  with  large  differences 
in  heat  flow,  by  convective  versus  conductive  means.  The  Russian 
group  of  Nazhukin  et  al.(1980)  dealt  with  stiffness  that  occured 
through  large  spacial  temperature  differences  created  through 
laser  irradiation  of  targets  and  interaction  with  the  resulting 
plasma  above  the  surface.  This  type  of  stiff  heat  problem 
commonly  comes  from  non- isothermal  chemical  kinetics  coupled  in  a 
model  with  transport;  refer  to  the  sections  in  this  review  on 
combustion  and  on  reactors. 

Mention  also  should  be  made  of  the  stiff  set  of  ODEs  that 
result  from  discretization  of  the  PDEs  that  describe  unsteady 
heat  transport.  See,  for  example  Bushard( 1976)  ,who  solves  the 
heat  conduction  equation  with  the  method  of  lines.  Wood (1977) 
discusses  the  solution  of  the  stiff  equations  that  result  from  a 
finite  element  discretization  of  the  heat  conduction  equation. 

Distillation  is  a  chemical  engineering  heat  transport 


process  that  can  result  in  stiff  numerical  models  because  of 
differences  in  the  liquid  hold-up  in  the  big  boiler  at  the  bottom 
of  the  column  compared  to  the  much  smaller  hold-up  on  the  plates. 
If  dynamics  of  the  vapor  traffic  is  included  together  with  that 
of  the  liquid,  stiffness  occurs  throught  the  great  difference  in 
heat  capacities  between  the  liquid  and  the  gas.  Tyreus  et 
al.(1975)  examined  stiffness  in  a  specific  model  of  a 
distillation  tower  and  found  it  became  more  severe  the  more 
difficult  the  separation  (from  high  purity  requirements  or  from 
the  components  to  be  separated  being  similar  in  their 
.volatilities).  An  adaptive  semi-implicit  Runge-Kutta  algorithm 
was  used  by  Prokopakis  and  Seider(1980)  in  a  model  in  which  the 
rapidly  changing  liquid  flow  rates  were  decoupled  in  a  sense  from 
the  relatively  slowly  changing  mole  fractions.  See  also 
Seider( 1982) .  Ozoe  et  al.  tackle  a  stiff  thermoacoustic 
convection  problem. 

This  type  of  stiffness  from  heat  capacity  differences  very 
commonly  occurs  in  reactors  which  contain  two  or  more  phases.  For 
example,  a  reactor  tube  containing  a  solid  particle  packing  that 
is  processing  a  gas  will  experience  stiffness  if  the  dynamics  of 
the  temperature  change  in  the  solid  and  gaseous  phases  are 
included  in  the  model. 

See  also  Churchill ( 1982)  for  a  review  of  stiff  heat 


transport  problems 


CHEMICAL  KINETICS 

This  is  by  far  the  largest  stiff  application  area.  Stiffness 
is  caused  in  the  vast  majority  of  cases  simply  by  a  large 
difference  among  the  reaction  rate  constants.  The  larger  the 
system  or  the  more  detailed  the  model,  the  more  likely  that 
stiffness  will  occur.  If  the  elementary  reactions  are  known,  the 
"law”  of  mass  action  dictates  the  form  of  the  rate  expressions: 
either  linear  in  a  concentration  variable  or  quadratic. 

Several  investigators  have  made  use  of  this  special  simple 
structure  of  mass  action  kinetics.  Edsbergt 1974) ;  Edsberq( 1976) ; 
and  Edsberg( 1982)  make  the  problem  set-up,  Jacobian  evaluation 
simplified  for  the  user  and  efficiently  handles  the  Jacobian,  but 
Dahlquist  et  al.(1980)  feel  more  can  be  done  to  make  use  of 
structure  as  well  as  the  users'  knowledge  of  the  stiffness.  This 
knowledge  often  consists  of:  a  partitioning  of  variables  into 
"stiff  variables"  and  "non-stiff  variables.";  a  fast  transient 
that  occurs  initially  only;  inherent  tight  stability  of  the  stiff 
variables  whose  concentrations  must  not  be  negative.  Karasalo  and 
Kurylo(1980)  point-out  an  advantage  in  keeping  these 
concentrations  artifically  non-negative  when  using  a  version  of 
GEAR.  Robertson ( 1975,  1976)  also  suggests  some  structure-related 
handling  of  the  Jacobian  for  faster  convergence. 

Packages  exist  specifically  for  the  mass  action  kinetics 
form.  Uhlen( 1979)  describes  KINRATE  and  KINBOX.  Edelson( 1976) ,  of 
Bell  Laboratory,  presents  a  simulation  language  and  compiler  for 
mass  action  kinetics;  he  uses  a  version  of  GEAR  in  forming  a 


package  called  BELLCHEM.  Rider(1977)  offers  CAKE,  user  friendly 
version  of  GEAR  that  makes  use  of  the  typical  sparse  structure  of 
large  mass-action  kinetic  equations.  Deuflhard  et  al.(1981)and 
Bader  et  al.(1982)  describe  LARKIN  to  handle  large  systems  of 
kinetic  equations.  Gottwald( 1981 )  gives  us  KISS  for  coupled 
chemical  reactions.  Stabler  and  Chesick( 1978)  have  also  written  a 
program  for  reaction  rate  equations  using  a  version  of 
GEAR. David ( 1977)  describes  a  FORMAC  program  for  direct 
integration  using  formula  manipulation  and  a  Taylor-made 
numerical  method;  Kennedy  and  Moore(1977)  also  recognized  the 
virture  in  using  a  Taylor-series  expansion  as  the  basis  for  a 
numerical  method  with  such  simple  functions. 

Enriqht  and  Hull (1976)  compare  numerical  methods  for  stiff 
kinetic  problems  and  found  that  the  backward  differentiation 
methods  were  superior  to  most  other  methods,  including  an 
implicit  Runga-Kutta  technique;  pitfalls  in  generalizing  such 
conclusions  are  explained  by  Enright( 1982 ) . 

The  steady-state  approximation  has  been  extensively  used  to 
eliminate  stiffness  in  chemical  kinetic  systems(Aiken  and 
Lapidus,  1975).  Noyes (1978)  discusses  the  importance  of  including 
reversible  reaction  when  the  approximation  is  made. 

Sensitivity  analysis  is  becoming  an  effective  means  of 
determining  appropriate  model  detail.  Koda  et  al.(1979)  studied 


automatic  sensitivity  analysis  of  kinetic  mechanisms  and 
developed  PAST  (Fourier  amplitude  sensitivity  test).  Dougherty 
and  Rabitz(1980)  look  at  the  sensitivity  of  hydrogen  combustion. 
Rabitz  gives  an  overview  of  this  area  applied  to  chemical 
kinetics.  Hwang(1982)  has  a  means  for  nonlinear  sensitivity 
analysis  in  chemical  kinetics.  See  also  Kuchel ( 1980) ,  Sundaresan 
and  Amundson ( 1980) ,  Dougherty  et  al.(1979),  Dove  and 
Raynor ( 1979) ,  Dickinson  and  Gelinas( 1976) ,  and  Lowe  et  al.(1977). 

A  specific  kinetic  application  area  not  covered  within  any 
of  the  other  Sections  of  this  review  is  pyrolysis.  Hautman  et 
al.(1981)  mention  that  at  low  conversions  the  primary  reactions 
govern  the  dynamics,  but  at  higher  conversion  the  secondary 
reactions  do.  Layokun  and  Slater(1979)  model  a  free  radical 
mechanism  of  propane  pyrolysis  and  solve  it  with  a  semi-implicit 
trapezoidal  rule.  A  number  of  thermal  cracking  models  were  solved 
in  detail  by  Sundaram  and  Froment( 1978) .  Liquid  phase  pyrolysis 
of  1,2  diphenylethane  was  studied  by  Miller  and  Stein(1981). 

An  overview  on  the  computational  techniques  for  the  study  of 
reaction  processes  is  avaiable  from  Edelson( 1981 ) .  A  stiff  model 
for  chemistry  in  interstellar  clouds  is  advanced  by  Prasad  and 
Huntress( 1980) .  Ross(1977)  mentions  the  problem  of  loss  of 
detailed  balance  when  applying  the  steady-state  approximation  and 
does  this  instead  by  a  Markov  matrix  method.  Rosenbaum  classifies 
certain  numerical  methods  as  "conservative"  (satisfying  the 


detailed  balance)  or  not.  Ong  and  Mason(1976)  discuss  a  different 
type  of  stiffness  in  kinetic  systems  than  the  type  we  have  been 
implicitly  focusing  on:  one  in  which  the  right  hand  side  is  the 
difference  of  two  large  terms;  they  convert  the  initial-value 
problem  into  a  two  point  boundary-value  problem  for  the  case  that 
the  right-hand  side  passes  through  zero.  See  also  the  entire 
volume  81,  number  25  pp  2309-2559  of  the  Journal  of  Physical 
Chemistry ( 1977) . 

LASERS 

Lasing  results  from  creating  a  highly  excited  vibrational 
state  in  a  group  of  molecules ("pumping”  to  an  inverted  energy 
state).  Then  a  remarkable  fact  of  nature  dictates  that  whatever 
the  "relaxations"  back  to  lower  energy  states  that  begin,  the 
resulting  photons  will  stimulate  other  of  the  still  excited 
energy  states  to  relax  in  the  same  way,  creating  an  autocatalytic 
effect  -  and  coherent  radiation.  These  relaxations  are  very  fast 
and  one  source  of  stiffness  in  the  modelling  of  lasers.  Shampine 
and  Gear (1979)  point  out  that  the  fast  pump-emission,  pump- 
emission  cycle  of  the  various  energy  levels  is  a  re-occuring 
stiffness  that  cannot  be  dealt  with  effectively  with  typically 
available  automatic  stiff  packages.  This  is  because  fast 
transients  occur  through-out  the  solution  so  accuracy  in  handling 
these  transients  can  be  continually  important.  Cukier  and 


Levine(1978)  mention  use  of  a  steady-state  approximation  during 
the  lasing  action,  but  details  on  this  apparently  ad  hoc  approach 
were  not  given.  There  is  an  initial  fast  transient  from  the  onset 
of  pumping  (Milonni,  1977),  that  is  at  t*0. 

Cukier  and  Levine  also  examine  the  sensitivity  of  the  full 
solution  for  a  model  of  an  HP  chemical  laser  and  find  only  a  few 
of  the  rate  constants  are  responsible  for  the  computed  gain.  A 
more  detailed  model  of  an  HP  laser  is  given  by  Ben-Shaul  and 
Peliks( 1979)and  by  Kerber  et  al.(1977).  It  should  also  be 
mentioned  here  that  the  rate  consants,  particularly  the  fastest 
ones  are  only  known  very  approximately,  errors  of  several  orders 
of  magnitude  are  not  uncommon. 

High  energy  lasers  usually  use  a  gas  as  the  lasing  medium 
for  best  efficiency  and  as  a  flowing  medium  to  remove  heat.  The 
addition  of  flow  can  require  the  laser  model  to  include 
hydrodynamics  coupled  to  the  chemistry.  The  time  scales  of 
transport  versus  reaction  are  much  different  here  (see  also 
sections  on  reactors  and  combustion).  Inclusion  of  hydrodynamics, 
translational  and  rotational  energy  interactions,  wall  effects 
and  the  like  are  necessary  to  actually  predict  the  perforamnce  of 
the  laser  -  to  compute  how  much  cower  output  can  be  extracted  and 
its  nature.  Additional  sources  of  stiffness  can  result  from  the 
variety  of  time  scales  amonst  the  three  very  different  types  of 


energy  transitions  for  a  molecule:  vibrational,  rotational, 
transulational  (see  the  section  in  this  review  on  molecular 
dynamics).  A  molecule  in  a  particular  energy  state  of  the  many 
potential  combinations  is  considered  a  distinct  species,  so  that 
a  large  number  of  highly  interactive  species  can  result  from  only 
a  few  different  molecules. 

A  meeting  conviened  at  the  U.S.Air  Force  Army  Weapons 
Laboratory  brought  in  stiff  experts  Shampine,  Gear,  Liniger, 
Hindmarsh,  and  Byrne  to  lecture  on  stiffness  and  hear  several 
talks  by  laser  modelers.  Some  of  these  have  been  referenced 
above,  others  are:  Franklin( 1977 )  who  spoke  on  modelling  general 
kinetic  processes  in  lasers;  Lundstrom( 1977) ,  of  the  Naval 
Weapons  Laboratory,  on  modelling  the  CO  laser;  Holmes(1977)  on 
C02/N2  vibrational  kinetic  equations;  Hines  on  CW  C02  electric 
discharge  modelling;  and  Young  and  Boris(1977),  of  the  Naval 
Research  Laboratory,  on  general  numerical  techniques  for  chemical 
kinetics  with  reactive  flow. 

Plasma  chemistry  has  been  studied  as  it  relates  to  laser 
discharge  and  target  interactions.  The  chemical  reactions  can 
include  neutral  molecule-neutral  molecule  collisions  as  well  as 
electron-ion,  electron-molecule,  and  molecule-ion  collisions. 
Roberts ( 1979)  presents  his  program  PLASKEM  for  this  problem. 
Pert(1978)  calculates  ionization  in  rapidly  changing  plasmas  in  a 
model  that  includes  hydrodynamics  but  integrates  these  two 


regimes  separately  in  a  steady-state  type  of  approximation.  The 
Russian  team  of  Mazhukin  et  al.(1980),  in  a  numerical 
investigation  of  laser  breakdown  of  a  dense  gas,  regect  a  version 
of  GEAR  and  the  method  of  lines  for  that  of  the  Russian 
Samarskii( 1971 ) ;  his  method  also  is  a  type  of  steady-state 
approximate  decoupling  of  kinetics  from  transport.  Christiansen 
and  Winsor(1980)  study  a  numerical  model  for  laser  targets, 
essential  to  the  feasibility  investigation  of  laser  fusion. 

Refer  to  Lawton  et  al.(1979)  for  numerical  work  on  the  high- 
pressure  infrared  xenon  laser;  Greene  and  Brau(1978)  for  RrP  and 
ArF  lasers;  Barker(1980)  for  infrared  multiphoton  decomposition; 
Pirkle  et  al.(1974)  for  pulsed  DP-C02  transfer  lasers;  and  to 
Bui (1979),  Bui (1980),  and  Bui (1981)  for  model  design  and  analysis 
of  a  new  type  of  blast-wave  induced  laser. 

MECHANICS 

The  term  "stiffness"  is  commonly  used  in  structural 
mechanics  in  a  much  different  sense  than  the  present  context.  The 
"stiffness"  matrix  results  from  a  linear  model 

f-Ky 

where  f  is  force  and  y  is  displacement.  The  denser  K,  the 
"stiffer"  the  problem  (finite  elements  used)  in  the  physical 
sense  that  there  are  more  interactions  among  the  components.  This 
would  imply  typically  a  more  stable  system,  but  a  stiff  system  in 


this  sense  is  not  necessarily  associated  with  fast 
transients. Large  second-order  systems  of  the  form 

My"  +  Cy’  +  Ky  «  f(t) 
y(0) ,  y'(0) ,  given 

occur  frequently  in  the  transient  analysis  of  dynamic  structures 
(Enright,  1980).  M,  C,  and  K  are  the  mass,  damping,  and  stiffness 
matrices,  respectively.  In  the  case  of  large  deflections,  the 
problem  becomes  nonlinear  through  K  depending  on  y.Here  the 
forcing  function  f(t)  can  make  accuracy  an  important 
consideration  for  any  component  of  y,  including  the  stiff  ones, 
at  any  time  in  the  transient;  the  numerical  stability  of  all 
components  is  essential.  This  special  type  of  stiffness  results 
in  computational  costs  very  dependent  on  requested  accuracy;  once 
a  step-size  has  been  selected  it  should  remain  relatively 
constant  throughout  the  solution,  there  being  no  boundary  layers. 
Thus  fixed  step-size  low-order  multistep  methods  have  been 
commonly  used,  although  Enright ( 1980)  has  pointed  out  the 
advantages  of  second-order  variable  step-size  approaches. 
Jensen( 1976 , 1974 )  of  the  Palo  Alto  Structural  Mechanics 
Laboratory  examines  stiffly  stable  third-order  methods  for  this 
same  problem.  See  also  Wright{ 1979 ) ,  Von  Pragenau( 1981 ) ,  and 
Addison( 1 980 )  for  the  linear  case  and  Park(1975)  for  the 
nonlinear  case.  This  last  work  is  particularly  interesting  as  it 
demonstrates  methods  unconditionally  stable  for  linear  problems 


are  not  so  for  the  nonlinear  case 


Jain  and  Jain(1981)  develop  hybrid  P-stable  methods  that 
improve  efficiency  for  solving  periodic  problems  in  celestial 
mechanics.  De  Silva  and  6rant{1978)  of  NASA  Ames  describe 
research  into  the  development  of  automatic  structrual  synthesis 
methods  for  turbine  disk  and  blade  assemblies.  Second  variation 
methods  resulted  in  systems  of  stiff  inhomogeneous  matrix 
Riccatti  equations. 

MOLECULAR  DYNAMICS 

Mathematical  models  of  atomic  and  molecular  dynamic 
interaction  have  successfully  predicted  macroscopic  physical 
properties  of  fluids.  Only  a  very  small  quantity  of  the  fluid  can 
be  modeled, to  limit  the  variables  to  hundreds  or  thousands  of 
particles.  Each  particle  can  theoretically  exert  a  force  on  all 
the  other  particles,  but  the  simplification  is  made  that  beyond 
some  cut-off  distance  the  force  is  too  small  to  consider  in  the 
model.  This  causes  a  discontinuity  in  the  interaction  potential 
that  is  a  source  of  error  in  many  molecular  dynamic  simulations. 

The  strength  of  the  interaction  is  a  very  strong  function  of 
distance.  Nearest  neighbors  are  thus  lead  by  a  rapidly  varying 
primary  force,  while  particles  farther  apart  change  more  slowly 
with  time.  This  natural  partitioning  by  distance  of  the  stiff 
from  non-stiff  variables  has  been  used  by  Streett  et  al.(1978) 


with  a  second-order  Taylor  series  method  and  extrapolation  to 
decrease  computation  time.  Their  application  was  to  108  methane 
molecules  with  five  sites  of  interaction  per  molecule  or  540 
potential  interactions.  Adequate  accuracy  was  established 
arbitrarily  on  the  basis  of  satisfaction  of  conservation  of 
energy  to  within  0.05  per  cent  per  1000  time  steps. 

The  mathematical  model  allows  "computer  experiments"  to  be 
performed  that  could  not  be  done  in  the  laboratory ,  or  would  be 
expensive  to  do.  Broughton  and  Abraham( 1980)  illustrate  this  in 
their  study  of  crystal-melt  interfaces.  They  use  a  variety  of 
GEAR. 

Heinzinger  et  al.(1978)  investigated  simulations  of  liquids 
with  ionic  interactions  and  found  higher  order  integration 
schemes  were  necessary  for  the  faster  rotational  motion  compared 
to  the  translational  motion.  Rossky  and  Karplus( 1979)  studied 
solvation,  Karplus  et  al.(1980)  studied  internal  dynamics  of 
proteins,  Rossky  et  al.(1979)  further  explored  solvent-solute 
interactions  all  using  a  version  of  GEAR  with  typical  time  steps 
on  the  order  of  E-16  seconds.  In  this  last  study,  the  limitation 
of  step  size  was  attributed  to  "rapid  liberational  motion  of  the 
water  molecules  and  the  correspondingly  rapid  change  in  the 
interaction  energy." 

Dove  and  Raynor { 1982 ) ,  Dove  and  Teitelbaum( 1979) ,  and  Dove 
and  Raynor{1979)  offer  an  interesting  approach  to  study  of 


vibrational  relaxation  in  hydrogen;  the  relaxation  was  treated  as 
a  chemical  kinetics  problem,  with  each  vibrational-rotational 
level  being  considered  a  distinct  species.  Reference  also  Gerlich 
et  a. (1980),  Haile  and  Graben( 1980) ,  and  Powles  et  al.(1979)  for 
other  stiff  molecular  dynamic  simulations. 

Microscopic  simulations  of  reacting  systems  has  evolved 
separately  from  molecular  or  atomic  dynamics  without  reaction. 
This  area  has,  in  addition  to  the  time  consuming  potential 
surface  evaluations,  multiple  times  scales  from  mechanical  versus 
the  chemical:  while  some  modest  time  frame  may  be  enou th  to  model 
motion,  a  reactive  interaction  is  typically  rare  on  that  time 
scale  -  and  fast.  See  Turner(1978)  for  a  review  of  reactive 
molecular  dynamics  and  consideration  of  how  the  interplay  of  the 
physical  and  chemical  on  the  molecular  level  can  influence 
macroscopic  physical  properties. 

NUCLEAR 

Safety  considerations  in  this  field  encourage  detailed 
dynamic  models  for  worst-case  numerical  experiments,  training, 
design,  and  control.  Thus  stiffness  can  be  identified  at  the 
atomic  stage,  mechanical  fuel-handling  stage,  spent-fuel  disposal 
stage,  or  the  overall  process  stage. 

The  radiolytic  decompostion  of  water  is  important  to  both 
moderator  and  coolant  chemistry  in  nuclear  reactors,  and  is  stiff 


because  of  disparity  of  rate  constants.  Carver  and  Boyd(1979)  and 
Boyd  et  al.(1980)  present  a  model  and  solve  it  by  their  own  stiff 
sparse  integration  method  for  mass  action  kinetics  called 
MACKSIM.  An  example  of  stiff  nuclear  reactions  are  given  in  Ward 
and  Fowler( 1980) ,  handled  by  a  GEAR  program. 

Lawrence  and  Doming ( 1977)  use  smoothing  and  extrapolation 
techniques  on  equations  stiff  because  they  include  the  greatly 
differing  time  constants  associated  with  prompt  and  delayed 
neutrons  in  subcritical  and  delayed  supercritical  transients  in  a 
fast  reactor.  The  so-called  "prompt  jump"  approximation  is  the 
same  as  the  commonly  used  steady-state  approximation;  Blenski  et 
al.(1978)  give  higher  order  singular  perturbations  for  this 
application.  Thermal  reactor  transients  are  moderately  stiff  and 
prompt  supercritical  reactivities  non-stiff. 

Devooght( 1980)  has  developed  a  more  general  steady-state 
method  for  nuclear  reactor  space-time  kinetics  as  used  for 
description  of  power  transients  in  fast  and  thermal  reactors.  He 
gives  a  review  of  stiffness  in  these  models  in  Devooght( 1982) . 

Rapid  ejection  of  a  control  rod  from  the  center  of  a  thermal 
reactor  is  solved  by  Carver  and  Baudouin( 1976)  using  a  version  of 
GEAR  with  their  own  sparse  matrix  solver.  They  also  solve  a  stiff 
test  model  of  a  control  rod  withdrawn  in  a  direction  parallel  to 
a  reactor  channel.  The  transient  is  a  long  200  seconds. 

Electronic  Associates,  Inc.,  of  West  Long  Branch,  N.J.,  are 


leaders  in  analog  computer  simulation  of  process  control  centers 
in  nuclear  power  plants  for  training  purposes.  The  real-time 
critical  nature  of  this  problem  and  the  curse  of  dimensionality 
for  digitial  simulation  explain  the  attractiveness  of  the  analog 
approach  here. 

Digital  nuclear  process  simulation  has  been  advanced  by  the 
Electric  Power  Research  Institute  (Bailey,  1982).  Patterson  and 
Rozsa(1980)  ,  of  Lawrence  Livermore  Laboratory,  describe  a 
nuclear  process  simulator  called  OYNSYL,  also  useable  for  more 
general  chemical  processes.  Chambers ( 1978)  documents  use  of  his 
AGR  (Advanced  Gas-cooled  Reactor)  digitial  simulator  for  real¬ 
time  solution.  Thompson  and  Tuttle(1982)  present  recent  software 
developed  at  Babcock  and  Wilcox  with  an  interesting  explanation 
of  its  historical  development  in  this  industrial 
environment. Halin( 1976)  discusses  performance  of  conventional 
stiff  methods  on  nuclear  problems  with  discontinuities. 

Borgonovi  et.  al(1980)  investigate  solution  of  stiff  models 
for  predicting  plutonium  inventory  on  a  continuous  basis. 
Carver(1981)  discusses  numerical  aspects  of  thermal-hydraulics. 
Gaf fney ( 1982) ,  of  Union  Carbide  Nuclear  Division,  surveys  methods 
for  solution  of  stiff  oscillatory  problems  as  arise  from 
magnetohydrodynamic  equations. 


PROCESS  INDUSTRIES 


Refer  to  the  reviews  of  Seider(1982)  and  Chen  and 
Schiesser( 1982) .  Also  to  Hylton(1982)  for  experience  with  CSNP. 
Barney  and  Johnson( 1975)  explain  the  incorporation  of  a  GEAR 
version  ordinary  differential  equation  solver  into  modular 
simulation  framework  of  DYNSYS; the  latest  version  of  this 
package,  developed  at  Lawrence  Livermore,  DYNSYL  is  described  in 
Patterson  and  Rozsa(1978).  Nilsen  and  Karplus( 1974 )  give  a  review 
of  continuous-system  simulation  languages.  Ockendon  (1980)  give  a 
survey  of  dynamic  simulation  of  Oxford  industry  problems. 

REACTORS 

Naturally  a  reactor  model  will  be  stiff  if  its  kinetic 
equations  are  stiff  or  if  there  is  a  difference  in  the 
characteristic  transport  time  from  the  reaction  time,  but  we  are 
concerned  here  with  other  types  of  stiffness  originating  from  the 
nature  of  the  reactor  model.  A  reactor  is  a  vessel  through  which 
reactant  continuously  flows  in  and  product  out;  there  may  be 
multiphase  flow  or  the  reactor  can  contain  a  fixed  solid  phase  on 
which  reaction  takes  place.  Introduction  of  each  new  phase 
increases  the  potential  for  stiffness,  depending  on  the  model 
detail,  as  different  phases  can  have  much  different  physical, 
chemical,  thermodynamic,  and  transport  properties.  For  example,  a 
tubular  reactor  containing  a  solid  particle  phase  changes  its 
temperature  much  more  slowly  than  a  gas  phase  flowing  through  it. 


The  strong  exponential  temperature  dependence  of  reaction  rates 
can  cause  very  sharp  temperature  spikes  as  one  marches  down  an 
exothermic  tubular  reactor.  This  nonlinearity  is  also  responsible 
for  the  existence  of  multiple  steady  states.  Aiken  and 
Lapidus( 1974)  give  an  example  of  a  non-isothermal  fluidized  bed 
that  is  very  stiff  and  posseses  three  possible  steady-states.  See 
also  Pan  et  al.(1979)  and  Michelson( 1976) ,  who  solves  this 
problem  with  a  semi-implcit  Runge-Kutta  method. 

Inclusion  of  a  diffusive  or  dispersion  term  under  conditions 
where  convection  very  much  dominates  the  flow  description  (high 
Peclet  numbers),  leads  to  stiff  computation  (see,  for  example, 
Shah  and  Parakos,  1975, and  Serth,  1975).  Smith(1980)  uses  a 
finite  element  approach  to  this  problem. Discontinuous  boundary 
conditions  for  tubular  reactors,  often  under  some  controversy, 
cause  stiffness. 

Interphase  mass  transport  with  reaction  can  produce  stiff 
boundary-value  problems  inherently  unstable  in  any  direction. 
Aiken(1982)  solves  a  gas  purification  model  of  simultaneous 
transport  of  two  gases  into  a  liquid  where  each  react  with  a 
third  species;  one  of  the  two  reactions  is  much  faster  than  the 
other. 

Karanth  and  Hughes(1974)  used  orthogonal  collocation  to 
solve  a  detailed  model  of  an  adiabatic  packed  bed  reactor, 
including  interphase  transfer  to  the  particles  and  intraphase 


particle  transfer.  Cavendish  and  Oh (1979)  solve  the  equations  for 
diffusion  and  reaction  in  a  bed  of  poisoned  automotive  catalysts 
pellets  by  applying  first  Galerkin's  method,  then  a  version  of 
GEAR.  Guertin  et  al.(1977)  use  exponential  collocation  on  some 
stiff  reactor  models .Rodrigues  and  Beira(1979)  and  Dias  et 
al.(1982)  model  and  solve  stiff  fixed-bed  adsorbers.  Eigenberger 
and  Butt(1976)  explain  a  technique  for  automatic  non-equdistant 
grid  size  space  for  finite  differences  on  reactors  with  steep 
gradients.  Edelson  and  Schryer( 1978)  compare  finite  difference 
with  finite  elements  for  one-dimensional  reactive  flow.  Varma  et 
al.(1976)  explore  a  number  of  computational  methods  for  tubular 
reactors.  Ramshaw( 1980)  discusses  the  use  of  the  steady-3tate 
approximation  in  reactive  three  dimensional  flow  .  Cho  and 
Joesph(1981)  solve  a  heterogeneous  model  for  moving-bed  coal 
gasification  reactors  and  remove  stiffness  with  the  steady-state 
approximation. Chin  and  Braun(1980)  solve  a  model  of  reacting  flow 
in  a  porous  medium;  George  and  Harris(1977)  of  in  situ  oil  shale 

retorting  -  all  quite  stiff  problems. 

See  also  the  section  in  this  review  on  combustion  where 

reactive  flow  in  the  confines  of  a  combustor  is  discussed. 
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INTRODUCTION 


The  user  would  like  to  obtain  the  solution  to  his  ordinary 
differential  equation  set  quickly  (in  both  his  time  and  in  the 
computer's  time)  and  in  a  convenient  form,  like  plots  or  a 
tabulation  with  automatic  choice  selection  of  where  or  when  to 
report  data. 

The  vast  majority  of  stiff  users  at  the  beginning  of  the 
1980's  use  early  version  of  Gear's  packaqe,  commonly  available 
through  users'  computer  centers.  By  numbers  of  users,  the  most 
common  problem  is  probably  small  (less  than  ten  or  twenty 
equations)  moderately  stiff  and  nonlinear.  It  only  needs  to  be 
solved  to  accuracy  of  about  one  per  cent  for  qualitative 
understanding  of  the  dynamic  relationships  among  the  variables. 

Necessary  information  from  the  user  obviously  includes 
functional  relationships  and  initial  conditions;  other  important 
but  possibly  expendable  instructions  would  be  accuracy 
requirements  and  integration  interval  (see  Klinzing( 1980)  for  a 
scheme  to  automatically  fiqure  out  this  interval).  Much  more 
additional  information  often  proves  to  increase  efficiency  of  the 
solution,  but  would  probably  not  be  welcome  to  the  average  user. 
This  includes  making  a  decision  on  whether  or  not  the  system  is 
stiff;  there  even  has  not  been  general  agreement  on  the  parts  of 
the  experts  on  what  stiffness  is  (Shampine,  1977; 1982).  In  this 
review,  stiffnes*  is  used  in  its  most  general  sense  of  "stubborn" 


numerical  method  behavior  where  differential  equations  are 
involved. 

The  user  interface  for  ODE  solvers  has  been  discussed  by 
Hindmarsh{ 1 978 ) ;  a  standard  was  developed  through  extensive 
discussions  with  individuals  from  six  DOE  laboratories.  This 
effort  was  intended  to  be  part  of  a  larger  effort  to  form  a 
standard  collection  of  ODE  solvers,  called  ODEPACK. 

Formation  of  an  ODE  package  as  automatic  as  possible  (but 
with  options  should  the  user  wish  to  specify) ,  winning  approval 
with  numerical  analysts  and  big  users,  then  the  public  is 
difficult  and  related  to  the  problems  associated  with  making 
comparisons  amonq  packaqes  (see  the  section  in  this  review  on 
this  topic).  See  Thompson  and  Tuttle( 1982)  for  a  description  of 
the  evolution  of  an  ODE  package  in  an  industrial  environment. 

Education  of  the  user  to  the  point  where  the  true  power  of 
an  approach  can  be  utilized  is  another  problem  area;  quite  often 
a  good  numerical  methods  course,  including  solution  of  ODEs  is 
not  in  an  engineers  curriculum.  See  "Getting  the  Power  to  the 
People"  by  Hindmarsh ( 1 980) ,  also  "What  Everyone  Solving 
Differential  Equations  Numerically  Should  Know,"  by  Shampine 
(  1980).  Also  Curtis(  1980) . 

Future  directions  to  go  in  the  development  of  general 
sol  >’are  is  discussed  in  Gear(  1982;  1981),  Cellier(  1982)  ,  and 

Dahlquist( 1980) .  This  includes  ultimate  automation  such  as 


automatic  method  selection  for  special  problems  like  stiff  and 
oscillatory  systems  dif ferential-algebraic  sets,  differential- 
difference  sets,  systems  with  discontinuities,  sparsity,  and 
such.  Better  interaction  of  experimentation,  differential 
equation  model  identification  and  parameter  estimation  with 
numerical  solution  requirements  needs  to  be  developed  (see  Koda 
et  al.,  1979,  for  an  automatic  sensitivity  analysis  program). 

METHOD  IMPLEMENTATION 

Packaqe  construction  from  a  method  of  interest  is  a  very  big 
task  with  a  number  of,  alas,  empirical  decisions.  Gear (1980) 
discusses  the  alchemy  side  of  software  development.  See  also 
sections  in  this  review  on  detecting  stiffness,  start-up,  on 
alqebraic  equations,  and  on  step-size/order  selection. 

Jackson  and  Sacks-Davis( 1980)  present  an  alternate 
implementation  of  variable  step-size  mutistep  formulas  for  stiff 
ODE.  Burrage  et  al.(1980)  discuss  implementation  of  singly- 
implicit  Runqa-Kutta  methods.  Ueberhuber( 1979)  suggests  an 
implementation  of  defect  correct  methods  for  stiff  ODE.  Sacks- 
Davis(1980)  implements  fixed  SD  formulas  with  fixed-leading 
coefficients  for  stiff  ODE.  Hindmarsh( 1 979 )  treats  software 
design  for  ODEs  from  PDEs. 


START-UP 


Gear(1980)  discusses  method  and  initial  step-size  selection 
in  multistep  ODE  solvers.  Early  automatic  codes  like  Gear's 
oriqinal  DIFSUB  required  the  user  specify  the  method (stiff  or 
non-stiff?  analytic  or  numerically-determined  Jacobian)  and  the 
initial  step-size. Usually  a  user  does  not  know  to  the  appropriate 
order  of  magnitude  what  size  the  first  step  should  be  and  may  not 
know  either  the  best  method  to  use.  The  easier  question  of  what 
method  is  the  more  important,  but  a  badly  guessed  initial  step 
can  waste  a  significant  amount  of  time  or  cause  one  to  miss  some 
interesting  behaviour.  Gear(1980)  examines  an  approach  that  uses 
the  initial  local  eigenvalues  to  make  the  (non-stiff)  method 
selection,  and  high  accuracy  estimates  at  nearby  points,  to 
select  a  good  initial  step  sequence  and  allow  a  high  order  to  be 
used  from  the  beginning  (where  it  can  be  particularly  valuable 
for  stiff  equations  and  problems  with  inital  zeros). 

Shampine( 1 978 )  studies  the  effect  finite  word  length  limits 
on  the  minimum  steo  size  has  on  solutions  containing  zeros, 
discontinuities,  and  stiff  transients.  He  finds  only  in.  this 
third  category  should  there  be  in  practice  limitations? 
suggestions  for  algorithmic  ways  out  are  given.  Thus  beginning 
with  extremely  small  steps  can  be  feasible  as  long  as  protective 
measures  are  present  to  avoid  over  zealous  increases  from  there 
skipping  interesting  phenomena.  Hindmarsh( 1977)  gives  reasons  why 


an  ODE  solver  may  take  steps  smaller  than  the  machine  precision 
of  representing  the  independent  variable. 

Gear(1980)  presents  Runge-Kutta  like  formulas  which  enable  a 
high-order  multi-step  method  to  be  used  from  the  start.  Only  one 
Runge-Kutta  step  is  needed  to  do  this  and  efficency  can  be  nicely 
increased  automatically.  The  Runge-Kutta  method  can  be  viewed  as 
an  exptrapolation  technique. 

STEP-SIZE/ORDER  CONTROL 

This  idea  enabled  Gear's  backward  difference  method  to  take 
a  "step  jump"  in  advancing  automatic  solution  of  stiff  (and  non¬ 
stiff)  ODE.  Gear  et  al.(1974)  and  Gear  and  Tu(1974)  examine  the 
effect  of  a  changing  mesh  size  on  the  stability  of  multi  step 
methods.  Skelboe( 1977)  studied  control  of  order  and  step-size  for 
multistep  methods  when  one  of  the  eigenvalues  is  close  to  the 
imaginary  axis;  an  instability  test  is  provided  to  automatically 
pick-up  on  when  this  is  the  case.  Lindberg( 1977)  characterizes 
the  optimal  step-size  sequence  for  stiff  methods.  See  also  the 
section  in  this  review  on  start-up. 

STIFFNESS  DETECTION 

The  automatic  detection  of  stiffness  is  related  to  the 

question  of  automatic  method  selection,  that  is,  the  choice  of  a 
stiff  or  non-stiff  method.  Gear(1980)  recommends  that  initially  a 


non-stiff  method  always  be  used.  Petzold( 1980)  and  Petzold( 1982) 
discusses  a  technique  that  uses  information  available  at  the  end 
of  each  step  to  make  a  decision  between  continuing  with  a  stiff 
or  a  non-stiff  method.  Shampine  and  Hiebert( 1977)  report  a  simple 
but  effective  method  for  following  occurance  of  stiffness  by  use 
of  the  Runga-Kutta  Fehlberg  (4,5)  formulas;  also  see 
Shampine( 1977) . 

Kennealy  and  Moore (1977)  show  a  heuristic  method  for 
detecting  stiffness  in  mass-action  kinetics.  Braekhus  and 
Aasen(1981)  explore  use  of  various  explicit  methods  for  detecting 
stiffness  in  problems  of  structural  mechanics.  See  too 
Gladwell ( 1 980 ) . 

Gordon  and  Shampine ( 1977)  mention  a  code  at  Sandia,  called 
DE,  that  solves  non-stiff  ODE  but  keeps  a  computer-eye  out  for 
the  occurence  of  stiffness.  These  authors  also  identify  two  other 
very  important  ways  of  detecting  stiffness,  non-automatically: 
based  on  physical  reasoning  and  based  on  computational  experience 
with  similar  equations.  There  is  another  non-automatic  way: 
through  proper  scaling  of  all  variables,  as  would  be  done  for  a 
singular  perturbation  analysis  or  in  preparation  for  putting  the 
equations  on  an  analog  computer.  See  Flaherty  and  0'Malley( 1979) 
for  an  automatic  scheme  for  this  scaling  on  a  digitial  computer. 

Shampine( 1980)  reports  a  special  definition  of  stiffness 
appropriate  for  implicit  A-stable  formulas;  this  definition  can 


be  recognized  using  information  already  available  during  the 
integration.  See  also  Sacks-Davis  and  Sharapine( 1981 ) . 

COMPARISON  OP  STIFF  METHODS 

Enright ( 1982)  explains  why  it  is  not  meaningful  in  general 
to  compare  different  methods  or  packages  containing  different 
methods  in  order  to  arrive  at  a  "best"  method.  However  such 
comparisons  can  point  to  weaknesses  in  a  method  or  code  (Enright 
et  al.,  1975;  Hull,  1980;  Enright  and  Hull,  1976).  As 
Enright  mentions,  it  makes  much  more  sense  to  compare  different 
implementations  of  a  given  method,  or  to  compare  two  related 
methods  or  packages.  Byrne  et  al.(1977)  compare  GEAR  and  EPISODE 
with  respect  to  appearance  to  the  user,  members  of  the  package, 
features  of  software  engineering,  and  the  basic  algorithms. 
EPISODE  performs  better  than  GEAR  for  problems  involving  waves  or 
re-occuring  stiffness,  but  GEAR  is  better  for  simple  decaying 
problems. 

Brown(1978)  offers  the  program  package  STIFF-DETEST  for 
comparison  of  stiff  ODEs.  See  Weimar  and  Clough(1979)  for  a 
critical  evaluation  of  the  semi-implicit  Runga-Kutta  methods  for 
stiff  systems.  Thompson( 1977)  and  Bushard( 1976)  have  performed 
comparisons.  See  also  Scherer ( 1976) ,  Carver  et  al.(1979),  and 
Chan  et  al. ( 1978) . 


TEST  PROBLEMS 


Packaging  considerations  make  it  hard  enough  to  compare 
methods ,  but  the  bottom  line  to  evaluation  of  software  is  how 
well  it  works  on  "typical"  stiff  problems.  Within  selected 
application  areas,  there  may  be  typical  non-linear  structure, 
degree  of  stiffness,  and  size,  but  in  general  there  is  not.  One 
reasonable  approach  is  to  examine  one  application  area  only,  or 
to  examine  many  different  test  problems,  representative  of  a 
spectrum  of  different  applications. 

whichever  approach  i*  taken,  there  is  virtue  in  consistency 
of  choice.  By  far  the  most  commonly  used  stiff  ODE  test  problem 
is  the  three  kinetics  rate  equation  set  of  Robertson ( 1975,  who 
cites  the  original  1966  article);  we  have  counted  at  least  23 
uses  of  this  equation  set  in  January, 1974-March,  1982.  Other 
favorites  are  given  in  Robertson( 1975) . 

Enright  et  al.(1975)  list  five  classes  of  stiff  problems 
with  a  number  of  examples  of  each  class:  linear  with  real 
eigenvalues,  linear  with  non-real  eigenvalues,  non-linear 
coupling (smooth  to  transient  and  transient  to  smooth),  non-linear 
with  real  eigenvalues  (most  mass  action  kinetics  problems  are 
here),  and  non-linear  with  non-real  eigenvalues;  eigenvalue 
ranges  are  given  for  most  of  these.  However  Shampine( 1977)  and 
Shampine  and  Hiebert ( 1977)  found  several  of  the  examples  did  not 
qualify  ,  by  their  definition,  as  being  stiff  -  although  some 


have  other  types  of  anoraolous  behavior.  Shampine( 1981 )  points  out 
a  number  oP  shortcomings  with  this  test  set  and  enumerates  ways 
of  improving  it. 

Enright  and  Hull (1976)  give  ten  test  problems  involving 
chemical  kinetics,  in  batch'  and  more  complex  reactors.  Johnson 
and  Barney(1976)  document  eleven  problems  they  used  for  testing. 
Hindmarsh  and  Byrne(1976)  give  a  diurnal  kinetics  problem  with 
the  re-occuring  stiffness  feature  characteristic  of  atmospheric 
reactions;  they  also  have  a  simple  diffusion-convection  problem 
for  use  with  the  method  of  lines.  A  variety  of  test  systems  are 
offered  by  Michelsen( 1976) ,  including  a  large  one  and  a 
differential-  algebraic  one.  Chan  et  al.(1978)  list  eight  stiff 
problems . 

There  are  also  available  several  testing  equations  of  a 
rather  special  nature.  Skelboe( 1977)  suggests  stiff  problems  of  a 
highly  oscillatory  nature,  as  does  Gaf fney( 1982) .  Fatunla( 1980) 
lists  six  examples,  some  of  which  are  stiff  and  highly 
oscillatory.  Oahlquist  et  al.(1980)  give  a  simple  stiff  equation 
with  a  turning  point,  a  stiff  nonlinear  oscillator,  and  a 
combustion  example.  Kreiss  and  Kreiss(1981)  consider  an  example 
of  a  stiff  two-point  boundary  value  problem.  Carver(1980)  has  two 
simple  hyperbolic  eauations  for  testing:  Burger's  equation  and  a 
model  for  a  counter-current  heat  exchanger. 


ALGEBRAIC  EQUATIONS 

There  has  been  little  comparative  advancement  during  recent 
years  in  solving  linear  or  nonlinear  equations  relevant  to  ODE 
solvers,  except  for  large  sparse  systems.  However  see  Hindmarsh 
et  al.(1978)  for  algorithmic  advancements  for  dense  linear  LU 
decomposition.  Johnson  and  Barney(1976)  test  five  conventional 
methods  for  solving  linear  algebraic  equations  (NINV,  SOLVE, 
DECOMP-SOLVE,  JACOBI,  AND  GAUSS-SEIDEL) .  Byrne(1976)  and  Byrne 
and  Hindmarsh( 1977)  consider  solution  of  linear  block  tridiagonal 
forms  arising  from  PDE  descretization.  Sherman  and 
Hindmarsh( 1980)  consider  solving  the  linear  equations  from  Newton 
iteration  on  a  nonlinear  sparse  set  by  the  YSMP,  Yale  Sparse 
Matrix  Package.  See  also  Hindmarsh ( 1977)  and  the  section  in  this 
review  on  sparse  systems. 

Shampine( 1979)  says  that  the  solution  of  the  algebraic 
equations  from  implicit  ODE  formulas  is  special.  He  found  that 
the  residual  was  the  appropriate  measure  for  acceptance  of  an 
approximate  solution;  a  way  to  do  this  and  the  advantages  are 
detailed. 

Hindmarsh ( 1977)  considers  the  idea  of  rank-one  updates  for 
the  inverse  of  the  Newton  iteration  matrix  in  the  context  of 
solvinq  stiff  ODEs,  but  the  results  are  disappointing. 
Enright( 1978)  gives  us  a  more  efficient  method  for  matrix 
factorization  after  a  change  in  step-size  or  order,  particularly 


qood  for  large  dense  systems.  See  Shampine( 1981 )  for  a  pertinent 
discussion  on  Jacobians  and  stiff  methods. Also  Eitelberg( 1979) . 

ANALOG  COMPUTATION 

There  are  some  advantages  for  use  of  an  analog  computer  for 

time-critical  computation,  that  is  for  some  real-time  needs  or 

for  stiff  computation.  There  has  not  been  a  widely-accepted 

recent  evaluation  of  this  usefullness,  however.  Reasons  why 

analogs  are  used  comparatively  infrequently  include  the  basic 
equipment  expense  (about  a  thousand  dollars  an  integrator)  to  buy 

an  analog  for  a  laboratory  that  typically  already  has  a  digital 

computer.  The  digitial  is  much  more  versatile.  Also  there  does 

not  yet  exist  "software"  in  the  digital  sense.  Therefore  one  must 

"patch"  an  analog  manually,  although  this  could  be  based  on  a 

diagram  written  by  someone  else.  Another  huge  inconvienence  is 

the  need  to  scale  the  problem  so  that  all  variables  vary  on  the 

same  normalized  interval  a  normalized  amount.  For  stiff  ODEs  this 

is  rather  equilvalent  to  requiring  the  user  set-up  his  problem  in 

dimensionlessized  singular  perturbation  form.  However  once  this 

is  done,  the  ease  of  parameter  variation  and  the  continuous 

graphical  availability  of  the  solution  on  an  oscilloscope  makes 

for  an  excellent  environment  to  throughly  explore  sensitivities 

of  parameters  and  interactions  of  variables  nonlinearly  related. 

Our  laboratory  owns  a  hybrid  EAI  PACER  1000/580  system  and 


an  BA12000  analog  (connected  serially  to  an  Apple  digial 
computer).  We  have  on-going  work  in  the  evaluation  of  the  analog 
for  stiff  computation,  as  occurs  from  mass-action  kinetic  models 

of  certain  combustion  reactions.  We  have  noted  an  upper  limit  on 
the  degree  of  stiffness  that  the  analog  can  handle  (for  example 
eigenvalue  spreads  greater  than  about  10X6  seem  impossible  to 
solve);  however  there  is  likely  to  be  a  relation  between  a 
variable  not  solvable  on  the  analog  and  a  one  that  is  really  not 
that  important  to  the  overall  solution.  But  this  is  often  evident 
once  the  ecruations  have  been  properly  scaled,  before  being 
patched. 

For  solutions  of  stiff  equations  on  a  hybrid  computer, 
reference  is  given  to  Kogan  et  al.(1980),  Karba  et  al.(1980), 
Neundorf ( 1981 ) ,  and  El-Zorkany( 1981 ) .  Stiff  problems  on  an  analog 
alone:  Bernard-Weil  et  al. ( 1978) .Refer  to  Gear(1977)  for  comments 
on  the  use  of  the  digital  for  real-time  dynamics. 

DIFFERENTIAL-ALGEBRAIC 

These  occur  in  models  of  power  systems,  control  systems  and 
from  application  of  the  steady-state  approximation  or 
perturbation  methods.  They  present  particular  problems  with 
determining  initial  conditions,  error  estimation  and  step  -  size 
selection  (Gear  et  al.1981).  If  the  algebraic  equation  resulted 
from  setting  a  derivative  to  zero,  singly  algebraic  equations  can 


sometimes  be  solved  explicitly  and  back-substituted,  for  example 
in  chemical  kinetics  of  the  mass  action  form.  If  this  elimination 
cannot  be  effected,  the  differential-algebraic  set  can  exhibit 
the  same  stiffness  as  the  original  fully-dif ferential  set.  Some 
problems  cannot  even  be  solved  with  stiff  methods,  without 
extensive  modifications  in  the  error  estimates  and  other 
strategies  of  the  code  ;  and  sometimes  they  apparently  cannot  be 
solved  at  all  by  stiff  methods (Petzold,  1981).  See  too  her 
package  DASSL  (Petzold,  1982). 

Liniger( 1979)  gives  us  multistep  and  one-leg  methods  for 
implicit  mixed  differential-algebraic  systems.  Soderlind( 1980) 
has  written  DASP3,  a  program  for  partitioned  stiff  ODEs  and 
differential-algebraic  sets.  Chua  and  Dew(1982)  attack  these 
mixed  systems  that  also  include  discontinuities.  Gross(1976) 
presents  a  method  that  makes  special  use  of  the  structure  in  the 
differential-algebraic  set; the  nonlinear  system  is  split  into  a 
stiff  part  with  a  sparse  Jacobian  and  a  nonstiff  part.  Datta  and 
Martens ( 1974 )  investigate  automatic  step  size  selection 
techniques  for  a  method  tailored  for  this  combination  of  equation 
types.  Refer  also  to  the  algorithm  of  Starner( 1976) . 


DIFFERENTIAL-DIFFERENCE 


Delay  terms  arise  in  lossless  transmission  line  modelling, 
in  ecomonic  modelling,  and  in  ecological  modelling,  to  name  only 
a  few  areas.  Brayton( 1974)  develops  conditions  for  numerical  A- 
stabilty  for  these  systems.  Bickart ( 1 981 )  offers  a  program 
package  for  differential-difference  systems.  Van  der  Staay(1982) 
explores  composite  intergrat ion-interpolation  methods. 
Bickart ( 1 982 )  determines  F-stable  and  F(alpha,  beta)-stable 
integration-interpolation  methods.  Weiderholt( 1976)  studies  the 
stabiltity  of  multi-step  methods  for  this  class  of  mixed  equation 
forms.  Carver(1977)  studied  the  efficient  handling  of 
discontinuities  and  time  delays  in  ordinary  differential 
equations.  See  also  Roth (1980),  Watanabe  and  Roth (1982),  and 
Moore ( 1 974 ) . 

DISCONTINUITIES 

Ellison( 1981 )  classifies  events  that  cause  discontinuities 
as  either  a  time  event  or  a  state  event.  Automatic  detection  of 
time  events  is  straight-forward,  detection  of  state  events  is  not 
but  is  achievable  on  examples  given;  an  integration  method 
schemes-up  with  the  detection  device  for  a  program.  Halin(1976) 
points  out  short  comings  of  popular  stiff  software  on 
discontinuities;  he  applies  a  "quasi-analytic"  integration 
technique.  De  Doncker( 1978)  presents  an  automatic  integration 


algoritm  (in  QUADPACK)  that  makes  use  of  a  nonlinear 
extrapolation  technique  to  jump  discontinuities.  Hay  et  al.(1974) 
also  have  a  means  for  detection  of  a  break  and  readjustment  of 
the  step  so  that  the  break  is  at  the  step's  end. 

GLOBAL  ERROR 

Users  often  do  not  realize  that  their  integration  package 
uses  their  requested  accuracy  requirement  to  match  against  ar 
estimate  of  the  local  error ,  not  the  actual  error  in  the  solution 
(global  error).  Lindberg ( 1977)  shows  for  stiff  problems  the 
advantage  of  keeping  the  global  error  at  the  maximun  allowable 
level  during  long  intervals.  Dahlquist( 1981 )  reports  work  in 
progress  to  extend  Lindberg' s  ideas  to  automatically  control  step 
size  on  the  basis  of  global  error;  application  is  made  to  a 
system  in  partitioned  form. 

Dew  and  West (1978)  consider  estimating  and  controlling 
global  error  in  Gear's  method.  Stetter( 1974)  considers  global 
error  estimation  for  non-stiff  problems;  Stetter( 1979 )  global 
estimation  in  Adams  predictor-corrector  codes. 


DECOUPLING 


In  addition  to  the  obvious  merits  of  smaller  size,  stiffness 
may  find  a  better  home  in  a  "decoupled"  or  semi-decoupled 
subsystem.  The  steady-state  approximation  is  the  best  known  way 
of  reducing  order,  but  the  differential-algebraic  set(see  this 
review)  may  be  just  as  stiff.  O'Malley  and  Anderson( 1979 )  discuss 
how  to  find  the  small  parameters  automatically  to  do  a  steady- 
state  approximation  (the  mathematically  sound  variety  obtained  by 
setting  a  parameter  to  zero  rather  than  a  derivative  to  zero) ; 
this  is  related  to  automatic  partitioning. 

Hofer  advocates  decoupling  stiff  from  non-stiff,  in  large 
systems  with  only  a  few  stiff  variables,  and  using  explicit 
techniques  on  the  non-stiff  part  and  implicit  methods  on  the 
stiff  part.  Enright  and  Kamel (1980)  study  selection  of  a  low- 
order  linear  model  using  the  dominant  mode  concept;  this  is 
related  to  lumping  and  modelling  questions. 

Nandakumar  and  Andres (1978)  explore  a  new  class  of 
algorithms  that  first  heuristicaly  decompose  large  systems  into 
groups  of  smaller  subsystems  that  share  similar  integration 
scales;  they  then  solve  individual  subsystems  and  combine 
iteratively.  Refer  also  to  the  decompositon  methods,  for  stiff 
equations,  of  Clasen  et  al.(1978),  Matthei j (  1 982 )  ,  Burka(  1982), 
and  the  problem-oriented  studies  of  Anderson( 1 980 )  for  control 
systems  and  Humpage  et  al. ( 1974)  for  power  systems. 


HIGHLY  OSCILLATORY  ODE 

Gaffney( 1982)  has  completed  a  critical  survey  and  testing  of 
software  (STRIDE,  BLEND,  STINT,  and  DIRK)  for  solving  stiff 
highly  oscillatory  ordinary  differential  equations;  none  of  these 
packages  are  given  very  high  marks  on  the  test  problem.  There  are 
also  such  things  as  highly  oscillatory  equations  that  do  not 
qualify  as  being  stiff  in  the  usual  sense  (no  large  negative  real 
eigenvalues),  although  they  can  be  "stubborn". 

Petzold( 1981 )  presents  a  numerical  method  for  this  (non¬ 
stiff)  highly  oscillatory  problem  as  does  Fatunla( 1980) .Gear  and 
Gallivant  1981 )  address  automatic  detection  of  highly  oscillatory 
behavior,  period  determination,  and  efficeint  integration.  See 
also  Gallivan( 1980)  and  Gear(1980). 

PARTITIONING 

The  practitioner  may  very  well  know  which  components  are 
highly  stable,  that  is  stiff,  and  which  are  not.  This  information 
can  be  used  to  make  the  numerical  solution  more  efficient  for  a 
variety  of  techniques.  Most  importantly,  such  knowledge  could  be 
used  to  make  a  modeling  simplification  to  remove  the  highly 
stable  component  from  the  model,  or  to  make  a  mathematical 
simplification  to  the  problem:  the  steady-state  approximation. 

This  approximation  has  the  terrific  property  of  being  better 
the  stiffer  the  system;  it  can  be  suprisingly  accurate  for  even 


weakly  stiff  systems.  It  is  tricky  to  apply  correctly  for  some 
systems,  however  (Aiken,  1982).  The  mathematical  basis  is  a  low- 
order  outer  approximation  in  singular  perturbation  theory  and 
this  can  be  quite  different  from  merely  mechanically  setting  a 
derivative  equal  to  zero.  A  few  recent  interesting  applications 
of  the  ad  hoc  version  of  the  approximation  are:  Chen  et 
al.(1979),  Devooght  and  Mund(1980),  Rao(1980),  Farrow  and 
Graedel ( 1 977 ) ,  Aronowitz  et  al.(1977),  Warner ( 1 977) ,  Cao  and 
Joesph( 1979 ) ,  and  sophisticated  versions  in  O'Malley  and 
Anderson( 1979 )  and  O'Malley  and  Flaherty( 1980) .  Application  to 
the  initial  conditions  can  eliminate  the  initial  transient  (Aiken 
and  Lapidus,  1975;  Alfeld,  1980). 

Soderlind( 1979)  discusses  some  stability  properties  of 
linear  multistep  compound  methods  on  a  system  partitioned  into 
two  sections.  Different  techniques  are  used  on  each  section. 
Palusinski  and  Wait(1978)  examine  methods  on  stiff  partitioned 
systems  into  one  linear  and  one  nonlinear  system  and  into  two 
nonlinear  systems.  Andrus(1979)  also  took  this  two  section,  two 
method  route.  See  also  Soderlind( 1980) ,  Soderlind( 1981 )  and 
Dahlquist( 1981 ) .  Enright  and  Kamel(1979),  Carver( 1982) ,  and 
Dahlquist  and  Fu-Fan(1982)  are  working  on  automatic  partitioning. 


SPARSE  SYSTEMS 


Large  systems  (more  than  ,  say  ,100  equations)  are  often 
sparse  because  there  are  usally  direct  interactions  among  only  a 
few  variables  in  the  set.  This  is  true,  for  example,  in  large 
kinetic  rate  equations.  Large  systems  resulting  from 
descretization  of  PDEs  are  sparse  with  special  structure.  Special 
handling  of  the  sparsity  can  both  reduce  storage  and  computation 
time.  Curt is (1977)  reviews  the  state  of  the  art. 

Because  the  efficency  of  Gear's  method  depends  heavily  on 
the  efficiency  of  matrix  operations,  particularly  for  larger 
systems,  Hindmarsh( 1 974 )  provided  a  more  versatile  package, 
called  GEAR  that  provided  several  matrix  options:  the  Chord 
method,  the  diagonal  method,  and  functional  iteration.  Later,  he 
added  an  option  for  banded  matrices,  as  from  PDEs,  in 
GEARS (Hindmarsh,  1975).  The  banded  structure  also  occurs  in 
models  of  stage-wise  processing  (Tyreus  et  al.,  1975).  Carver  and 
Baudouin( 1976)  used  this  package  to  solve  a  stiff  set  of  242  ODEs 
modelling  neutron  kinetics  and  transport;  they  found  that  only 
the  chord  method  allowed  the  solution  to  be  reached  in  reasonable 
time  but  storage  was  near  to  machine  capacity  and  20  seconds  were 
required  for  each  Jacobian  evaluation  and  decomposition.  They 
therefore  added  to  the  package  a  method  from  the  Harwell 
subroutine  library  for  solution  of  large  linear  equations,  which 
stores  only  non-zero  entries  and  uses  a  pivotal  scheme  optimal  in 


some  sense.  The  resultinq  package  is  called  FORSIM.  The  greatest 
savings  was  probably  in  the  numerical  approximation  to  the 
Jacaobian:  the  matrix  is  evaluated  first  to  find  non-zero 
elements;  then  the  Jacobian  is  evaluated  by  perturbing  as  many  of 
the  individual  variables  as  effect  only  one  derivative. 

Sherman  and  Hindmarsh( 1980)  present  GEARS ,  a  package  for 
stiff  sparse  ODE ,  using  YSMP  (Yale  Sparse  Matrix  Package).  The 
two  main  sparse  techniques  here  are  a  special  method  for 
computing  finite  difference  approximations  to  the  Jacobian  and 
YSMP  non-pivoting  Gausian  elimination  linear  equation  solver. 

Schaumberg  et  al.(1930)  and  Zlatev  et  al.(1980)  analyze 
implementation  of  a  Gustavson  storage  scheme  and  a  generalized 
Markowitz  pivotal  strategy  for  large  stiff  linear  ODEs. 
Enright ( 1979 )  also  examines  stiff  sparse  linear  ODEs,  and 
suggests  modifications  to  GEAR  for  three  classes  of  linear 
equtions  and  four  levels  of  structure. Guy  Rabbat  et  al.(1979) 
mention  sparse  matrix  techniques  have  allowed  time  domain 
analysis  of  circuits  with  hundreds  of  elements,  but  large  scale 
integrated  circuits  present  the  challenge  of  solving  thousands  of 
active  devices.  Johnson  and  Barney(1976)  look  at  several  sparse 
techniques  (SIMULT,  IMP,  and  LINEQ4)  . 

For  strategies  solving  applications  resulting  in  PDEs-turned 
ODEs,  see  Iserles( 1981 ) ,  Sincovec  and  Madsen( 1975) ,  Melgaard  and 


Sincovec( 1981 ) ,  Hunding( 1980)  and  Karasalo  and  Kurylo{ 1981 ) .  The 
many  other  example  implementations  of  sparse  techniques  include: 
Franke( 1980) ,  transient  field  problems;  Watson(1976)  for  CSNP 
III;  Gross(1976)  for  power  systems;  Carver  et  al.(1979)  for  mass- 
action  kinetics;  Dove  and  Raynor(1979,  1982),  molecular  dynamics; 
Enright ( 1980) ,  structural  mechanics;  Sincovec  et  al.(1981)  for 
describer  systems;  Prasad  and  Huntress( 1980) ,  interstellar 
clouds;  Atkinson  et  al.(1980)  ,  atmospheric  pollution;  and 
Thompson  and  Tuttle(1982)  for  process  problems. 

UNSTABLE  PROBLEMS 

Lindberg( 1974)  discusses  the  fact  that  many  stiff  methods 
fail  to  detect  inherent  instabilty  of  an  equation,  particularly 
when  larqe  negative  eigenvalues  turn  positive.  Aiken(1982)  notes 
automatic  methods  can  skip  over  an  explosion  in  the  model. 
Aiken(1982)  notes  a  very  common  and  very  "stubborn"  numerical 
problem  in  studying  selectivity  in  qas  purification  operations, 
related  to  the  occurence  of  positive  eigenvalues.  Hoppensteadt  et 
al.(1981)  propose  a  numerical  method  that  focuses  on  the  positive 
eigenvalues.  Brown(1978)  examines  the  error  behavior  of  multistep 
methods  applied  to  unstable  differential  equations. 


PACKAGES 


Most  recent  stiff  packages  generally  available  are  mentioned 
in  various  sections  of  this  review.  Table  I  presents  a  summary  of 
most  of  these.  Lawrence  Livermore  National  Laboratory  (LLNL)  has 
led  the  way  for  developing  general  and  special  purpose  stiff 
packages  for  the  user.  Outside  of  LLNL,  GEAR  can  be  obtained  from 
the  National  Energy  Software  Center(NESC) ,  Argonne  National 
Laboratory,  9700  South  Cass  Avenue,  Argonne,  Illinois  60439; 
identify  GEAR  as  NESC  No.  592. 


i 


NAME 


General 

DGRUNG 

EPISODE 

FACSIMILE 

GEAR 

GRK4T 

IMP 

LSODE 

SDBASIC 

STIFF  3 

STINT 

STRIDE 

TRAPEX 


TABLE  I.  STIPF  SOFTWARE 

+ - f - 

COMMENT 

+ - + - 


Two-stage  semi-implicit 
Runge-Kutta 
Re-occuring  stiffness 

Derived  from  DIFSUB 
Rosenbrock  methods 
Implicit  midpoint 
More  ease,flexiblity 
Second-derivative 
Semi-implicit  R-K 
Cyclic 
Implicit  R-K 
Extrapolation 
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Boness( 1979) 

Byrne( 1981 ) 

Curtis( 1980) 

Hindmarsh( 1974 ) 

Kaps  and  Rentrop( 1979) 
Lindberg( 1974) 
Hindmarsh( 1980) 
Enright( 1975) 
Michelsen( 1976) 

Tendler  et  al.(1978) 
Butcher,  et  al.(1979) 
Enright  et  al.(1975) 


Direct  Variants  of  GEAR 
DSTPGT 

GEARS  Sparsity 

GEARBI  2-3  dimensional  PDE 

GEARV  For  parallel  processors 


Thompson  and  Tuttle(1982) 
Hindmarsh( 1979) 

Hindmarsh( 1979) 

Morris  et  al.(1979) 


Special  Applications 

CAKE 

Kinetics 

CSDT 

PDE,  adjust  mesh 

DISPL2 

Collocation-B-splines 

FORSIM  VI 

Sparse  techniques 

K INRATE 

Kinetics 

KISS 

Kinetics 

LARKIN 

Large  kinetics 

SETKIN 

Kinetics  preprocessor 

Diferential-Algebraic 

DASP3 

DASSL 

Also  partitioning 

EPISODEIB 

Banded  Jacobian 

FAST 

Transulator 

GEARIB 

Banded  Jacobian 

GEMS 

Extension  of  IMP 

LSODEI 

Linearly  implicit 

Process  Simulation 

CSMP 

CSMP  III 

DPS 

HYNSYL 

Ridler ( 1977) 

Janac( 1978) 

Byrne( 1981 ) 

Carver( 1979) 

Edsberg( 1974) 

Gottwald( 1981 ) 

Bader  et  al.(1982) 

Dickinson  and  Gelinas( 1976) 


Soderlind( 1980) 
Petzold( 1982) 
Hindmarsh( 1979) 
Stutzman  et  al.(1976) 
Hindmarsh( 1979) 
Babcock  et  al.(1981) 
Hindmarsh( 1980) 


Hylton( 1982 ) 

Watson  and  Gourlay( 1976) 
Sebastian  et  al.(1981) 
Patterson  and  Rozsa(1978) 
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Subject  listing  in  alphabetical  order  (note  various  stability 
types  are  listed  by  heading  under  STABILITY. 

AVERAGING 

Miranker( 1982) ?  Hoppensteadt ( 1 979 ) ;  Miranker( 1978 ) ;  Persek 
and  Hoppensteadt ( 1978) ;  Hoppensteadt  and  Mir anker ( 1977 ) ;  Coleman 
et  al . ( 1977 ) ;  Liniger( 1976) . 

BLENDED 

Cash (1981,  1982);  Skeel  and  Kong (1977). 

BLOCK 

Bond  and  Cash (1979). 

COLLOCATION 

Burka(1982);  Finlayson( 1 982 ) ;  Dias  et  al.(1982);  Chin  et 
al .  ( 1 979 ) ;  Guertin  et  al.(1977);  Varma  et  al.(1976);  Wouk(  1976); 
Bushard( 1976) ;  Shah  and  Paraskos ( 1975) ;  Serth(1975);  Carey  and 

Pinlayson( 1975) ;  Michelsen  and  Villadsen( 1974) ;  Karanth  and 
Hughes( 1974) ;  Scholander  and  Svedberg ( 1974) ;  Brunner( 1974) . 

CONTRACTIVE 

Liniger( 1982) ;  Soderlind( 1981 ) ;  Mingyou( 1981 ) ;  Ruehli  et 

al.(1980);  Odeh  and  Liniger( 1980) ;  Dahlquist  and  Jeltsch( 1979) ; 
Nevanlinna  and  Liniger( 1979 ) ;  Dahlquist( 1979) ;  Nevanlinna  and 
Liniger{ 1978) . 


COMPOSITE  MULTI-STEP 

Tendler  et  al .(1978);  Lucey(1975);  8ickart( 1974) . 

CYCLIC 

Gupta( 1979);  Tendler  et  al.(1978);  Mihelcic( 1978) ; 

Michelcic( 1977) ?  Cash(1977). 

DEFECT  CORRECTION 

Ueberhuber( 1979) ;  Frank  and  Ueberhuber( 1977) . 

DELAY 

Bickart( 1982) ;  Cryer(1974). 

DISCONTINUOUS 

Chua  and  Dew(1982);  Struwe( 1981 ) ;  Tuttle( 1981 ) ;  O'Malley  and 
Flaherty( 1980) ;  Halin(1979);  Srivastava  et  al.(1979);  De 

Doncker{ 1978) ;  Mannshart ( 1 978) ;  De  Hooq  and  Weiss(1977); 
Carver( 1977) ;  Halin(1976);  Lake  et  al.(1976);  Luke  et  al.(1975); 

Macccomack  and  Paullay{ 1974 ) ;  Hay  et  al.(1974). 

EXPLICIT 

Fatunla( 1980 ) ;  Alf eld( 1 979 ) ;  Kubicek  and  Visnak( 1974) . 
EXPONENTIAL  FITTING 

Cash( 1981 ) ;  Iserles ( 1 98 1 ) ;  Cash(1981);  Raptist 1980) ; 

Rentrop( 1980) ;  DeGroen  and  Hamker( 1979) ;  Miller ( 1 979 ) ; 
Iserles( 1979 ) ;  Strehel  and  Peper(1979);  Miranker( 1978) ; 


Norsett( 1978) ;  Iserlesf 1978) ;  Rosenbaum( 1978 ) ;  Iserles( 1977) ; 
Murphy( 1977) ;  Gear(1977);  Jackson ( 1976 ) ;  Ehle(1975);  Sarkany  and 

Liniger( 1974) ;  Jackson  and  Kenue(1974);  Chu  and  Berraan( 1974) ; 
Makela  et  al.(1977);  Chun-Yat( 1974) ;  Meister ( 1974 ) . 

EXTRAPOLATION 

Duff  and  Nowak(1982);  Hoppensteadt  and  Miranker( 1979) ; 
Lawrence  and  Doming  ( 1977) ;  Gladwell(  1976) ;  Cash (1976); 
Lindberg( 1974) . 

FUNCTIONAL  FITING 

Iserles(  1977) . 

HYBRID  METHODS 

Fatunla( 1982) ;  Jain  and  Jain(1981). 

IMPLICIT  EULER 

Frank  and  Ueberhuber( 1977) . 


LINEAR  MULTISTEP 

Nolting( 1982) ;  Liniger ( 198 1 ) ;  Van  Veldhuizen( 1981 ) ; 
Butcher( 1981 ) ;  Fatunla( 1980) ;  Sinha(1980);  Gear(1980);  Warming 
and  8eam(1979);  Soederling( 1979) ;  Dahlguist( 1979) ;  Alfeld( 1979) ; 
Grigorieff  and  Schroll ( 1978 ) ;  Dahlquist( 1978) ;  Jeltsch( 1 978 ) ; 
Varah(1978);  Nevanlinna( 1 977) ;  DeHoog  and  Weiss(1977); 
Jeltsch( 1977 ) ;  Kong(1977);  Prothero( 1 976 ) ;  Liniger( 1975) ;  Gupta 


and  Wallace* 1975) ;  Williams  and  De  Hoog*1974). 

MOLTIDERIVATIVE 

Cash(1981);  Burrage* 1980) ;  Jeltsch* 1979) ;  Jeltsch* 1977) ; 
Brown(1977);  Jeltsch* 1 976 ) ;  Fuchs(1976);  Gennin* 1 974 ) . 

NONLINEAR  MOLTISTEP 

Lee  and  Preiser* 1978 ) ;  Lee* 1974). 

ONE  STEP 

Van  Brokhoven* 1980) ;  Cash* 1978);  Mannshart* 1978) ; 

Watanabe* 1978) ;  Cash* 1975);  Prothero  and  Robinson* 1974 ) ;  Van 
Veldhuizen* 1 974 ) ;  Brunner ( 1974 ) ;  Norsett ( 1 974 ) . 

ONE  LEG 

Liniger* 1982) ;  Sand* 1981);  Dahlquist  et  al.(1980); 
Brown* 1979);  Nevanlinna  and  Liniger* 1979) ;  Dahlquist* 1979) ; 
Liniger* 1979) ;  Dahlquist* 1 978) . 

ORDER  STARS 

Hairer ( 1979 ) ;  Wanner  et  al.(1978). 

OSCILLATORY  EQUATIONS 

Miranker ( 1 982 ) ;  Fatunla* 1 982 ) ;  Jain  and  Jain* 1981); 
Petzold* 1981 ) ;  Cash(1981);  Gourlay* 1980) ;  Ruehhli  et  al.(1980); 

Brown* 1980);  Fatunla* 1980) ;  Rinzel  and  Miller ( 1980) ;  Warming  and 


Beam( 1979 ) ;  Hoppensteadt ( 1 979 ) ;  Hoppensteadt  and  Miranker( 1979) ; 
Bui ( 1979 ) ;  Kreiss ( 1 979 ) ;  Auslander  and  Miranker( 1979) ;  Jain  and 
Jain( 1979 ) ;  Ruehli  et  al.(1978);  Miranker  and  Veldhuizen( 1978 ) ; 
Miranker( 1978) ;  Jeltsch( 1978) ;  Hoppensteadt  and  Miranker( 1977) ; 
Skelboe( 1977) ;  Smith (1977);  Amdursky  and  Ziv(1977);  Lambert  and 
Watson( 1976) ;  Miranker  and  Wahba(1976);  Fatunla( 1 976) ; 
Gupta(1976);  Snider  and  Fleming( 1974) . 

QUADRATURE 

Iserles( 1981 ) . 


RUNGE-KUTTA 


Implicit 

Cash(  1982) ;  Zlatev( 1 98TT;  Hufndsdorfer  and  Spijker(  1981 ) ; 
Mingyou( 1981 ) ;  Cash(1981);  Burrage( 1979) ;  Van  der  Houwen( 1979 ) ; 
Bui  and  Bui(1979);  Burrage  and  Butcher ( 1979 ) ;  Crouzeix( 1979) ; 
Scherer ( 1979 ) ;  Dahlquist  and  Jetsch( 1979) ;  Eitelberg ( 1979 ) ; 
Curtis( 1979) ;  Butcher ( 1 979 ) ;  Varah(1979);  Freidli ( 1978) ; 
Burrage( 1978) ;  Palusinski( 1978) ;  Iserles ( 1 978 ) ;  Alexander* 1977) ; 
Bickart( 1977) ;  Butcher ( 1976) ;  Fuchs(1976);  Cash(1975); 
Ehle(1975);  Ehle  and  Lawson( 1 975 ) . 


Cooper ( 1979 ) . 


Semi-explicit 


Semi-implicit 

Cash( 1982) ;  Prokopakis  and  Seider ( 1981 ) ;  Weimer  and 
Boness{ 1979) ;  Kaps(1979);  Clough ( 1 979 ) ;  Cash(1979);  Bui(1979); 


Kaps  and  Rentrop( 1 979 ) ;  Burrage( 1978 ) ;  Lapidus  et  al.(1978); 
Freidli( 1978) ;  Cash(1976);  Michelsen( 1976) . 

SINGLY  IMPLICIT 

Burrage( 1980) . 

SECOND  DERIVATIVE 

Sacks-Davis  and  Shampine( 1981 ) ;  Sacks-Davis( 1980) : 

Enright{ 1978) ;  Gupta( 1978) ;  Sacks-Davis ( 1 977 ) ;  Kennealy  and 
Moore(1977);  Hill(1976);  Brown(1976);  Kubicek  and  Visnak( 1974) ; 

Enright( 1974) . 


SECOND  ORDER  EQUATIONS 


Addison( 1980) ; 

Odeh 

and 

Liniger( 1980) ; 

Heinrich  and 

Zienkiewicz( 1979) ; 

Van 

der 

Houwen( 1979)  ; 

Hairer(  1979) ; 

Gear(1978);  Jensen( 1 976 ) ;  Liniger  and  Gagnebin( 1974 ) . 

SEPARABLY  STIFF 

Lambert ( 1981 ) . 

SINGULAR  PERTURBATION 

Matthei j ( 1982) ;  Mattheij  and  0'Malley( 1982) ;  Kreiss  and 
Kreiss( 1981 ) ;  Petzold( 1981 ) ;  Soerderlind  and  Dahlquist( 1981 ) ; 
Smith( 1981 ) ;Sanchez-Palencia  and  Lobo-Hidalgo( 1 980 ) ;  DeVooght  and 
Dahlquist  et  al.(1980);  Mund(1980)j  Kahil  and  Kokotovic( 1980) ; 

Barton( 1980) ;  O'Malley  and  Flaherty( 1 980 ) ;  Brandt( 1979) ; 
Auslander  and  Miranker( 1979 ) ;  De  Groen  and  Hamker( 1979) ;  Heinrich 


and  Zienkiewicz( 1973) ;  Michell  and  Christie( 1979) ;  Hsiao  and 
Jordan( 1979 ) ;  Matthei j ( 1979) ;  Miller( 1979) ;  Flaherty  and 
O'Malleyt 1979) ;  Bourgeat  and  Tapiero( 1 979 ) ;  Andrus( 1979) ; 
Hoppenstesdt  and  Miranker( 1 979) ;  Come(1979);  Hoppensteadt ( 1979) ; 

Kreiss( 1979 ) ;  Miranker( 1978) ;  Emel'yanov( 1978) ;  Persek  and 
Hoppensteadt ( 1978) ;  David(1977);  Flaherty  and  0'Malley( 1977) ; 

Robertson ( 1975) ;  Aiken  and  Lapidus( 1975) ?  Flaherty  and 

0'Malley( 1975) ;  Dontchev( 1974) ;  Kreiss( 1974) ;  Aiken  and 
Lapidus( 1974) . 

SPLINES 

Rentrop( 1980) ;  Hill ( 1976) . 

STABILITY 

General 

Butcher ( 1981 ) ;  Lambert ( 1 980 ) ;  Brown(1979)j  Jury(1978); 
Dahlquist{ 1978) ;  Glaser ( 1 978 ) ;  Bickart  and  Jury(1978); 

Jury ( 1977) ;Dahlquist( 1976 ) . 

A-stable 

Bui  and  PoonH981);  Zlatev(  1981 ) ;  Iserles(  1981 ) ;  Odeh  and 
Liniger( 1980) ;  Tadeusiewicz  ( 1980) ;  Van  Brokhoven( 1980) ? 

Wanner( 1980) ;  Galantai( 1980) ;  Burrage  and  Butcher( 1979) ; 
Bui(1979);  Warming  and  3eam(1979);  Cooper  and  Sayfy(1979); 
Scraton( 1979) ?  Kaps(1979);  Tripathy  and  Rao(1978);  Bickart  and 

Jurv(1978);  Dahlquist ( 1 978 ) ;  Iserles ( 1 978 ) ;  Jeltsch ( 1 978 ) ;  Wanner 


et  al.(1978);  Watanabe( 1978 ) ;  Lee  and  Preiser( 1978) ; 
Triqiante( 1 977) ;  Brown(1977);  Butcher ( 1 977 ) ;  Jeltsch( 1977) ; 
Lombardi ( 1977) ;  Fuchs(1976);  Liniger ( 1976) ;  Chipman( 1976) ; 
Wanner ( 1976) ;  Jackson( 1976) ;  Jeltsch( 1976) ;  Cash(1976); 
Cash(1975);  Ehle  and  Lawson{ 1975) ;  Butcher( 1975) ;  Ehle(1975); 
Prothero  and  Robinson( 1974 ) ;  Williams  and  De  Hoog(1974); 
Marzulli( 1974) ;  Brandon( 1 974 ) ;  8rayton( 1974) ;  Liniger  and 
Gaqnebin( 1974) ;  Norsett ( 1974) ;  Axelsson( 1974) ;  Genin(1974); 
Nevanlinna  and  Sipila( 1974 ) . 


A( alpha)  Stable 

Galantai( 1980) ;  Kaps{1979);  Bickart  and  Jury(1978); 

Grigorieff  and  Schroll( 1978) ;  Michelcic( 1978) ;  Jeltsch( 1977) ; 


Liniger( 1975) . 


A(alpha,r)  Stable 

Nolting( 1982) . 

A ( 0 )  Stable 

Rodabaugh  and  Thompson( 1979 ) ;  Freidli  and  Jeltsch{ 1978) ; 
Jeltsch{ 1976) ;  Liniger(  1975)  . 

AO  Stable 

Jeltsch( 1976) . 

An  Stable 

Zlatev( 1981 ) . 

Algebraically  Stable 

Burrage{  1978) . 


Mihelcic( 1977) 


Almost  A  Stable 


Krogh( 1981 ) 


Asymptotically  Stable 


B  Stable 

Hundsdorfer ( 1981 ) ;  Burrage  and  Butcher { 1 979 ) ; 

Crouzeix( 1979 ) ;  Scherer( 1979) ;  Jeltsch( 1979) ;  Dahlquist  and 
Jeltsch( 1979) . 

D  Stable 

Veldhuizen( 1981 ) . 

P  Stable 

Bickart( 1982) . 


G  Stable 

Dahlquist( 1978) y  Nevanlinna( 1976) . 

I  Stable 

Wanner  et  al.(  1978).~* 

L  Stable 

Scraton( 1981 ) ;  Day?1980);  Cash(1980);  Lambert ( 1980) ; 

Bui (1979);  Scraton ( 1 979 ) ;  Bui(1979);  Eitelberg ( 1 979 ) ; 

Fatunla( 1978) ;  Bui{1977);  Trigiante( 1977) . 

L(alpha)  Stable 

Cash{ 1980 ) . 

L ( 0 )  Stable 

Gourlayt 1980) . 

Ln  Stable 

Zlatev( 1981 ) . 


Nonlinear  Stability 
Dahlquist( 1982) ;  Soederlind ( 1 98 1 ) y 

Butcher( 1980 ) ;  Burrage( 1980) ;  Brown(1979); 

Butcher( 1979 ) ;  Dahlquist( 1978) ;  Burrage(1978 

Whiworth( 1 978 ) ;  Rodabaugh  and  Thompson ( 1978 ) ; 


Burrage  and 
Burrage  and 
;  Cooper  and 
Trigiante( 1977) ; 


Liniger ( 1977) ;  Nevanlinna( 1977) ;  Wanner( 1976) ;  Liniger  and 
Odeh( 1976) ;  Dalquist( 1975) ;  Butcher( 1975) ;  Lambert ( 1974) ; 
Brandon ( 1974 ) . 

P  Stable 

Bickart( 1982) ;  Patunla( 1982) ;  Jain  and  Jain(1981); 
Cash( 1981 ) ;  Jain  and  Jain(1979). 

S  Stable 

Day(1980);  Alexander ( 1 $77 i ;  Verwer ( 1977) ;  Prothero  and 
Robinson( 1974 ) . 

Stiffly  A  Stable 

Ehle  and  Lawson ( 1975) . 

Stiffly  Stable 

Nolting (1982);  Watkins ( 1981 ) ;  Jeltsch( 1979) ; 

Watanabe( 1978) ;  Albrecht ( 1978) ;  Jain  and  Srivastava( 1978 ) ; 
Tendler  et  al.(1978);  Varah(1978);  Jeltsch( 1977) ;  Rao  and 

Iyengar ( 1976) ;  Jeltsch( 1976) ;  Gupta  and  Wallace( 1975) ;  Prothero 
and  Robinson( 1974) ;  Jensen( 1974) j  Bickart  and  Rubin(1974). 

Strong  A  Stable 

Watanabe( 1978) . 

Strongly  Stable 

Struwe(  1981 ) ;  Taubert ( 1 976 )  ?  ^Gear  and  Waanabe(  1974 ) ; 
Lee( 1974) . 

Strong  Stiffly  Stable 
Watanabe( 1978)1 

Zero  Stability 


Zlatev( 1978 ) 


TURNING  POINTS 

Ponisch  and  Schwetlick( 1981 ) ;  Moore  and  Spence( 1980) ; 
Miranker  and  Morreeuw( 1974 ) . 

TWO  STEP 

Iserles( 1981 ) ;  Dahlquist  et  al.(1980);  Odeh  and 
Liniger( 1980) ; 

UNSTABLE  EQUATIONS 

Aiken( 1982) ;  Hoppensteadt  et  al.(1981);  Mazurkin( 1980) ; 
Kreiss( 1979) ;  Serth(1975);  Lindberg( 1974) . 
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